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Abstract 


The  time-dependent  behavior  of  a  transversely  magnetized,  two- 
dimensional  plasma:wall  sheath  has  been  studied  through  particle 
simulations,  with  the  aim  of  modelling  plasma  behavior  in  the  vicinity 
of  the  limiters  and  walls  of  magnetized  plasma  devices.  The  simula¬ 
tions  have  shown  that  the  cross-field  sheath  between  a  wall  and  a 
plasma  is  a  turbulent  boundary  layer,  with  strong  potential  fluctua¬ 
tions  and  anomalous  particle  transport.  The  driving  mechanism  for 
this  turbulence  is  the  Kelvin- Helmholtz  instability,  which  arises  from 
the  sheared  particle  drifts  created  near  the  wall.  Provided  it  is  re¬ 
plenished  by  an  internal  flux  of  particles,  the  sheath  maintains  itself 
in  a  dynamic  equilibrium,  in  which  the  linear  edge  instability,  the 
nonlinear  dynamics  of  the  particles  and  the  outward  particle  diffusion 
all  balance  each  other.  The  sheath^assumes  an  equilibrium  thick¬ 
ness  of  order  lx  ~  5  p,  .  ancPmaintains  large,  long-lived  vortices,  with 
amplitudes  6<t>  ~  2T,/e, ‘which  drift  parallel  to  the  wall  at  roughly 
half  the  ion  thermal  velocity.  The  sheath  also  maintains  a  large, 
spatially-averaged  potential  drop  from  the  wall  to  the  plasma,  with 
A <t>  »  -l.5T,/e~,  in  sharp  distinction  with  the  unmagnetized  sheath, 
where  the  plasma  potential  is  higher  than  at  the  wall.  Accompanying 
the  long- wavelength  vortices  is  a  spectrum  of  shorter-wavelength  fluc¬ 
tuations,  which  extend  to  |k|  p,  ~  1  and  u  ~  uci,  and'which  induce  an 
anomalous  cross-field  transport,  A  central  result  is  that  the  anoma¬ 
lous  transport  scales  like  Bohm  diffusion,  at  least  when  u>p,  >  2u>ct.  At 
lower  densities,  u >p,  <  2u;cl,  the  diffusion  coefficient  has  an  additional 
factor,  proportional  to  the  density.  These  results  enable  us  to  model 
the  cross-field  sheath  by  a  simple  boundary  condition,  which  relates 
the  particle  flux  through  the  sheath  to  the  edge  density  and  which 
can  be  used  as  input  in  any  model  designed  to  obtain  the  bulk  plasma 
properties. 


Contents 

1  Introduction 

2  Previous  Work 

3  The  Simulation  Model 

4  Transient  Behavior  of  the  Sheath 

4  1  Evolution  of  the  Total  Number  of  Particles  .... 

4.2  Evolution  of  the  Electrostatic  Potential  . 

4.3  Evolution  of  the  Fourier  Modes . 

4.4  Minimum  Length  for  a  Two- Vortex  Steady-State  . 

4.5  Saturation  Mechanisms:  Climax  and  Coalescence  . 

4.6  Behavior  when  k\\  ^  0 . 

5  Fluid  Theory 

5.1  The  Nonlinear  Cross-Field  Equations . 

5.2  Growth- B  '  °s  of  the  Kelvin-Helmholtz  Instability  .  . 

6  Steady-State  Structure  of  the  Sheath 

6.1  Simulation  Results . 

6.2  Dependence  of  the  Sheath  Thickness  on  System  Size 

6.3  Localized  Solutions  for  the  Steady-State  Vortex  . 

6.4  Periodic  Solution  for  the  Steady-State  Vortex  .  .  .  . 

6.5  Power  Spectra . 

6.6  Source  of  the  Short- Wavelength  Fluctuations  .  .  . 

7  Transport 

7.1  Orbits  of  Test  Particles  . 

7.2  Dynamics  without  Electron  Gyromotion . 

7.3  Numerical  Estimates  of  the  Diffusion  Coefficient 

7.4  Estimate  of  Discrete-Particle  Effects . 

7.5  Analytic  Estimate  of  the  Diffusion  Coefficient  .  .  . 

7.6  Scaling  of  the  Diffusion  Coefficient  with  Density 


8  Summary 


9  Conclusions 
Acknowledgments 

A  The  One-Dimensional  Sheath 

A.l  Analytic  Expressions  for  the  Electric  Field . 

A. 2  Comparison  with  one-dimensional  particle  simulations 

References 


•V-  A m  .\A  vV  J»V -Av-  ««  *  A  ^  A  .VAj  V ^  r_i  "lA-a  **  ’A* ^jIAaTla 


wo 


,**  4 

A 


.  yv. -r«  %  ^ 


1  Introduction 

The  following  is  a  report  on  the  results  of  our  particle  simulations  of  the 
magnetized  plasma-wall  sheath.  This  is  a  study  of  plasma  transport  per¬ 
pendicular  to  a  magnetic  field,  in  a  plasma  bounded  by  a  conducting  wall. 
The  objective  is  to  model  plasma  behavior  in  the  vicinity  of  limiters  and 
walls  of  magnetized  plasma  devices.  Our  approach  has  been  to  use  our 
two-dimensional,  bounded  particle  simulation  code  ES2[1],  as  a  tool  for 
the  investigation  of  edge  effects,  in  an  idealized  model  which  retains  the 
essential  features  of  the  edge  plasma. 

Our  simulations  have  shown  that  the  cross-field  sheath  between  a  wall 
and  a  plasma  is  not  a  static  structure,  but  is  in  fact  a  turbulent  boundary 
layer,  with  strong  potential  fluctuations  and  anomalous  particle  transport. 
The  driving  mechanism  for  this  turbulence  is  the  Kelvin-Helmholtz  insta¬ 
bility  which  arises  from  the  sheared  particle  drifts  created  near  the  wall  by 
the  strongly  non-neutral  sheath.  Provided  it  is  replenished  by  an  internal 
flux  of  particles,  coming,  for  instance,  from  a  central  bulk  plasma  or  from 
a  diffuse  ionization  of  neutrals,  the  sheath  will  maintain  itself  in  a  dynamic 
equilibrium,  in  which  the  linear  edge  instability,  the  nonlinear  dynamics  of 
the  particles  and  the  outward  particle  diffusion  all  balance  each  other.  It 
is  important  to  emphasize  that  the  turbulent  behavior  of  the  sheath  is  a 
completely  spontaneous  phenomenon,  which  arises  from  the  self-consistent 
plasma-wall  interaction,  and  which  does  not  require  the  imposition  of  ex¬ 
ternal  fields  for  its  sustenance. 

We  have  found  that  the  cross-field  sheath  assumes  an  equilibrium  thick¬ 
ness  of  order  lx  ~  5  p,,  and  that  it  maintains  large,  long-lived  vortices,  with 
amplitudes  S(f>  ~  2Ti/e,  which  drift  parallel  to  the  wall  at  roughly  half  the 
ion  thermal  velocity.  The  sheath  also  maintains  a  large,  spatially-averaged 
potential  drop  from  the  waill  to  the  plasma,  with  A<£  «  —1.5 Ti/e,  in  sharp 
distinction  with  the  unmagnetized  sheath,  where  the  plasma  potential  is 
higher  than  at  the  wall.  Accompanying  the  long- wavelength  vortices  is 
a  continuous  spectrum  of  shorter-wavelength  fluctuations,  which  extend  to 
|k|  pi  ~  1  and  u>  ~  u and  which  induce  am  anomalous  cross-field  transport 
of  particles.  A  central  result  of  this  paper  is  that  the  anomalous  transport 
in  the  sheath  scales  very  much  like  Bohm  diffusion,  at  least  when  cjp,  >  2u>ct. 
At  lower  densities,  such  that  u>p,  <  2wCi,  the  diffusion  coefficient  is  found  to 
have  an  additional  factor,  proportioned  to  the  density.  These  results  enable 
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us  to  model  the  entire  cross-field  sheath  by  a  simple  boundary  condition, 
which  relates  the  particle  flux  through  the  sheath  to  the  edge  density.  Tins 
boundary  condition,  which  measures  the  sheath  impedance  to  particle  flow, 
is  usable  for  any  model  designed  to  obtain  the  bulk  plasma  properties,  and 
in  which  the  detailed  sheath  dynamics  are  unimportant. 

As  an  aside,  we  would  like  to  briefly  describe  our  simulations  in  more 
general  terms,  by  using  concepts  of  non-equilibrium  thermodynamics.  The 
vortices,  which  behave  as  convection  cells,  can  be  considered  dissipative 
structure$[2],  into  which  the  plasma  has  rearranged  itself  so  as  to  maximize 
heat  and  particle  fluxes  to  the  boundaries.  The  vortices  have  formed,  of 
course,  subject  to  the  constraints  inherent  in  the  physical  system;  these 
are  energy  and  momentum  conservation,  and  the  physical  length  and  time 
scales  available  to  the  plasma.  The  source  of  free  energy  winch  maintains 
these  structures  resides  in  the  temperature  difference  between  the  reservoir 
of  hot  ions  and  electrons,  which  are  continuously  fed  into  the  system,  and 
the  wall,  which  is  kept  at  absolute  zero,  insofar  as  it  is  a  perfect  absorber 
of  incoming  particles.  Our  system  is  analogous  to  the  configuration  lead¬ 
ing  to  Rayleigh- Benard  convection [3},  where  convection  cells  and  enhanced 
transport  are  also  driven  by  a  temperature  gradient.  However,  the  detailed 
dynamics  of  each  system,  and  the  physics  of  the  media  considered,  are,  of 
course,  very  different. 

The  work  described  in  this  Memorandum  has  also  been  presented  in 
[4],  A  videotape  of  a  computer  animation  showing  some  of  the  dynamic 
features  of  the  cross-field  sheath[5],  including  the  time  evolution  of  potential 
surfaces  and  density  profiles,  can  be  obtained  from  the  Plasma  Theory  and 
Simulation  Group. 


2  Previous  Work 


The  properties  of  the  cross-field  sheath  can  be  understood  as  arising  from 
the  combination  of  two  physical  effects,  the  formation  of  a  static  plasma- 
wall  sheath,  and  the  existence  of  a  sheared-flow  instability.  These  effects 
have  previously  been  studied  in  completely  distinct  physical  situations.  As 
both  fields  of  study  are  of  relevance  to  the  present  work,  we  shall  briefly 
review  the  work  which  has  been  done  in  both  of  these  hitherto  separate 
domains. 

Let  us  first  consider  static  models  of  the  sheath.  The  original  work  on 
unmagnetized  sheaths  goes  back  to  Tonks  and  Langmuir [6],  but  the  study 
of  magnetized  sheaths  is  more  recent.  Particularly  close  to  the  situation 
considered  in  the  present  paper  are  the  works  of  Delbege  and  Bein[7]  and 
Chodura[8],  in  which  the  authors  have  considered  the  magnetized  sheath 
which  forms  between  a  wall  and  a  plasma.  However,  there  are  crucial  differ¬ 
ences  between  their  work  and  the  content  of  the  present  paper:  their  anal¬ 
ysis  is  one-dimensional,  and  furthermore  assumes  a  magnetic  field  which 
is  oblique,  and  not  strictly  parallel,  to  the  wall.  In  this  configuration,  the 
electron  dynamics  are  dominant,  and  the  magnetized  sheath  is  indeed  qual¬ 
itatively  similar  to  the  unmagnetized  sheath. 

Let  us  now  consider  the  second  field  of  study,  that  of  shear-flow  instabil¬ 
ities  in  plasmas.  There  is  a  large  literature  going  back  to  the  original  work 
of  Buneman[9]  on  the  “Diocotron”  instability.  Somewhat  more  recent  work 
includes,  for  instance,  that  which  is  presented  in  articles  of  Gould[lO]  and 
Knauer[ll].  The  plasma  Kelvin-Helmholtz  instability,  as  opposed  to  the 
one-species  Diocotron  instability,  has  been  treated  by  Chandrasekhar [12] 
as  an  MHD  instability.  In  the  present  paper  however,  we  have  especially 
drawn  on  work  by  Byers[l3],  Miura  and  Pritchett[14],  and  Pritchett  and 
Coroniti[l5].  We  have  also  noted  the  work  of  Horton  et  al.[16],  in  which 
the  Kelvin-Helmholtz  instability  is  studied  in  a  configuration  with  similar¬ 
ities  to  the  present  one.  We  shall  discuss  these  similarities,  as  well  as  basic 
differences  between  their  model  and  ours,  in  the  summary  of  Section  8. 


3  The  Simulation  Model 


The  physical  configuration  embodied  in  our  computer  simulation  is  an  ideal¬ 
ized  model  of  the  tangential  edge  conditions  which  might  exist  in  a  tokamak 
or  in  a  variety  of  plasma  magnetrons.  We  assume  that  the  essential  plasma 
dynamics  are  two-dimensional,  with  all  motion  confined  to  the  plane  per¬ 
pendicular  to  the  magnetic  field,  with  magnetic  field  lines  exactly  tangential 
to  the  boundary.  The  magnetic  field  is  homogeneous  and  shear-free,  and  the 
model  is  restricted  to  electrostatic,  collisionless  modes.  These  assumptions 
imply  that  in  all  that  follows,  we  shall  be  concentrating  on  the  behavior  of 
flute-like  electrostatic  modes,  with  k\\  =  0,  driven  by  charge  separation,  and 
localized  in  the  edge  layer  In  a  real  device,  we  are  focussing  on  a  distance 
much  smaller  than  the  scale-length  of  any  magnetic  shear 

By  focussing  on  flute-like  modes,  we  are  insuring  that  the  assumption 
of  a  magnetic  field  strictly  parallel  to  the  wall  does  not  correspond  to 
a  singular  configuration.  Rather,  we  are  assuming  that  the  system  self- 
consistently  selects  an  instability  with  sss  0,  for  which  the  assumption  of 
two-dimensional  geometry  is  automatically  valid. 

We  modelled  the  two-dimensional  plasma-wall  system  outlined  above 
with  an  electrostatic  particle-in-cell  simulation  code,  the  computer  program 
ES2[1].  This  two-dimensional  program  is  of  the  explicit  type,  in  that  it 
incorporates  full  electron  and  ion  dynamics,  and  simulates  plasma  behavior 
on  the  time-scale  of  the  electron  gyro-period,  u;"1  [17]. 

The  simulation  region  is  shown  in  Fig.(l).  The  particles  move  in  the  (x- 
y)  plane,  in  which  they  are  subject  to  the  external,  perpendicular  magnetic 
field  Bo,  and  in  which  they  also  interact  with  each  other,  according  to 
their  self-consistent  electrostatic  fields.  Note  that  both  particle  velocities 
and  coordinates  are  restricted  to  two  components,  with  v  =  (Ux,t>y)  and 
x  =  (x,y). 

The  simulation  region  is  rectangular  in  shape,  with  periodic  boundary 
conditions  on  the  particles  and  fields  in  the  direction  parallel  to  the  wall 
(y),  with  the  wall  defining  the  left-hand  boundary.  The  wall  is  taken  to  be  a 
perfect  conductor,  imposing  a  constant  potential  at  x  =  0,  and  is  modelled 
as  a  perfect  absorber  of  impacting  electrons  and  ions,  with  no  reflection 
or  reemission  of  particles.  When  a  particle  hits  the  wall,  it  is  lost  to  the 
inner  plasma,  and  its  charge  is  immediatelv  accounted  for  as  distributed 
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surface  charge,  which  in  turn  is  a  boundary  source  in  Poisson’s  equation. 
In  solving  Poisson’s  equation,  the  program  also  automatically  accounts  for 
all  charges  induced  on  the  wall  by  the  particles  moving  inside  the  plasma. 
In  our  simulations,  the  wall  is  kept  “floating”;  it  is  not  connected  to  an 
external  circuit,  and  charge  flow  can  occur  only  from  the  plasma  to  the 
wall. 

The  boundary  conditions  at  the  right-hand  side  of  Fig.(l)  are  an  at¬ 
tempt  to  simulate  as  best  as  possible  a  semi-infinite  plasma  in  j*  >  Lx. 
The  boundary  conditions  are  the  following.  For  the  fields,  the  surface  at 
x  =  Lx  is  an  equipotential.  For  the  particles,  an  inversion  symmetry  condi¬ 
tion  is  imposed[18];  a  particle  exiting  at  y  with  velocity  ( vx,vy )  is  returned 
at  y'  —  Ly  —  y  with  inverted  velocity  (~vx,  -vy).  Thus,  the  boundary  at 
x  —  Lx  can  be  thought  of  as  a  very  fine  wire  mesh  which  shorts-out  the 
tangential  electric  fields,  but  which  is  transparent  to  particles.  The  mesh 
separates  two  plasma  regions,  which  are  exactly  symmetrical  through  the 
point  (Ir,£y/2),  and  which  consist  of  the  actual  simulation  plasma,  and  a 
virtual,  spatially  inverted  “twin”  plasma.  This  inversion-symmetry  condi¬ 
tion  differs  somewhat  from  the  one  outlined  in  [18],  because  the  boundary 
at  x  —  Lx  is  taken  to  be  an  equipotential.  We  believe  that  this  condition 
shields  the  actual  simulation  region  from  the  fields  in  its  virtual  twin,  and 
thus  reduces  the  interference  between  the  two  regions. 

In  ES2,  we  have  also  modelled  a  distributed  plasma  source,  with  which 
ions  and  electrons  are  continuously  created  in  the  system.  This  is  done  by 
creating  electron-ion  pairs  spatially  at  random,  and  at  a  constant  temporal 
rate  (a  constant  number  of  pairs  is  created  at  each  time  step).  The  electron 
and  ion  in  each  pair  are  initially  created  on  top  of  each  other,  and  are  given 
random  and  uncorrelated  velocities,  chosen  from  Maxwellian  distributions 
with  Tf  —  Tt.  This  source  of  plasma  can  be  thought  of  as  resulting  from 
the  ionization  of  a  uniform  background  of  neutral  atoms,  and  is  essential  in 
maintaining  the  system  in  a  steady-state  on  the  time-scale  of  the  cross-field 
transport.  Over  shorter  time-scales,  an  initial  loading  of  warm  electrons 
and  ions  with  no  further  pair  creation,  is  sufficient  to  produce  the  turbulent 
cross-field  sheath,  with  a  slow  net  diffusion  of  particles  to  the  wall. 

In  Fig.(l).  we  have  suggested  the  presence  of  an  underlying  heat  en¬ 
gine,  which  drives  the  processes  in  the  edge  layer.  We  have  done  tins  by 
indicating  teir  ltures  for  the  wall  and  for  the  plasma  source:  the  wall, 
which  absorbs  ail  mcoi  particles  and  emits  none,  can  be  thought  of 
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as  a  heat  sink  at  absolute  zero;  the  plasma  source  on  the  other  hand,  is 
an  infinite  reservoir,  maintained  at  a  constant  temperature  T  /  0.  In  our 
fictitious  “heat  engine”,  it  is  this  temperature  differential  between  wall  and 
the  plasma  source  which  provides  the  free  energy  necessary  for  maintaining 
turbulent  behavior  in  the  edge  layer. 


Run  1:  Parameters  for  the  Initial  Conditions. 


mi/mt 

40 

Vtc 

i 

Ax 

1 

s' 

1.005 

Vti 

0.158 

Ay 

1 

u>ce 

1 

Pc 

1 

At 

1 

Wpe 

0.2 

Pi 

6.32 

Nx 

33 

Up, 

0.032 

Lz/p, 

5.06 

256 

0.025 

40.5 

Ucxtmaz  ~  1500 

Upi/ujCi 

1.15 

Cl 

ft 

II 

o. 

5 

Nt(t  = 

0)  15000 

0.182 

■ai 

56 

N,(t  = 

0)  15000 

Table  1:  Parameters  for  Run  1. 


4  Transient  Behavior  of  the  Sheath 


We  shall  now  consider  the  results  of  an  ES2  particle  simulation,  with  pa¬ 
rameters  given  in  Table  1.  This  computer  run,  which  we  shall  refer  to  as 
Run  1,  is  the  central  example  to  be  discussed  in  this  paper.  We  initialized 
the  system  with  a  plasma  uniform  up  to  the  wall,  and  ran  the  simulation 
with  a  constant  creation  of  particles.  After  a  fairly  lenghty  transient  phase, 
the  system  settled  into  a  quasi-steady-state,  in  which  the  outward  transport 
of  particles  balanced  the  influx  from  the  source,  and  in  which  the  qualita¬ 
tive  structure  of  the  electrostatic  fields  remained  the  same.  In  the  present 
and  following  sections,  we  shall  discuss  the  transient  phase  of  evolution 
of  the  plasma-wall  sheath,  and  examine  the  growth  and  saturation  of  the 
electrostatic  fields.  In  particular,  we  shall  present  analytic  results  for  the 
linear  growth  of  the  Kelvin-Helmholtz  instability,  and  we  shall  qualitatively 
describe  the  mechanisms  underlying  the  subsequent  nonlinear  evolution  of 
the  system.  A  discussion  of  the  steady-state,  as  characterized  by  its  fields, 
spectra,  and  particle  transport,  is  deferred  to  Section  5  and  to  the  following 
sections. 


Let  us  examine  the  conditions  implied  by  the  parameters  displayed  in 
Table  1.  In  Run  1,  the  plasma  region  was  initially  filled  with  N  =  15000 
particles  of  each  species,  uniformly  distributed  in  space  and  initialized  with 
Maxwellian  velocity  distributions  with  equal  temperatures,  T,  =  Te.  The 
ion-electron  mass  ratio  used  in  the  simulation  was  m,/me  =  40.  With  the 
normalizations  assumed  in  ES2  (e0  =  1,  me  =  e),  and  with  e  =  me  — 
0.018  and  BQ  =  1,  we  have  in  the  units  of  ES2  uce  =  1,  uct  =  0.025, 
vte  =  ( Te/me yt2  =  1,  and  vti  =  (T./m,)1^2  =  0.158.  The  initial  particle 
densities  were  ne  =  n,  =  N/LxLy  =  1.83,  such  that  u)pe  =  0.182  and 
u ip,  =  0.0287,  giving  the  dimensionless  ratios  uipe/uJce  =  0.182  and  u>p,/u>c,  = 
1.15.  We  also  have  the  values  for  the  gyroradii  and  the  Debye  lengths, 
pe  =  1,  p,  =  6.32,  and  A^e  =  A^,  =  5.5,  so  that  the  physical  lengths  satisfy 
Pe  •C  A je  =  Xj,  <  pi.  The  dimensions  of  individual  cells  on  the  mesh 
were  Ax  =  Ay  =  1,  so  that  we  could  resolve  wavelengths  down  to,  but 
not  including,  the  electron  gyroradius.  Our  choice  of  the  initial  number 
of  particles  was  such  as  to  ensure  a  large  number  of  particles  per  Debye 
square,  with  n\2de ,  =  55.  This  large  value  of  the  plasma  parameter  was 
needed  to  keep  the  thermal  fluctuations  of  the  system  small  compared  to 
the  amplitudes  of  the  collective  modes. 

In  Run  1,  our  choice  of  a  fairly  low  particle  density,  with  uip, / u/c,  =  1.15, 
enabled  us  to  run  ES2  for  very  long  times  at  a  reasonable  computational 
cost,  and  thus  to  accumulate  a  large  amount  of  simulation  data.  In  Section 
7.6  we  shall  discuss  the  results  from  simulations  with  higher  densities,  but 
shorter  running  times.  However,  Run  1  exhibits  all  the  basic  features  of 
the  cross-field  sheath,  and  we  use  it  as  our  central  example. 

The  overall  system  size  was  Lx  =  32  and  Ly  =  256,  with  a  grid  .Yx  x 
Ny  =  33  x  256.  Thus  the  system  was  elongated  parallel  to  the  wall,  with 
Ly/ Pi  =  40.5  and  Lx/pt  =  5.06.  Our  choice  of  Ly  >  Lz  was  dictated  by  our 
desire  to  accomodate  several  unstable  wavelengths  of  the  initially  dominant 
Kelvin-Helmholtz  mode. 

The  simulation  was  run  with  a  constant  background  creation  of  electron- 
ion  pairs,  at  an  average  rate  of  =  1.005  pairs  injected  into  the  system 
per  time  step  (At  =  1).  We  can  better  appreciate  the  injection  rate 
if  we  define  defines  a  characteristic  “filling-up"  time,  tj,u  =  X/stn]:  this 
is  the  time  in  which,  without  outward  loss,  the  number  of  particles  in  the 
system  would  double.  Thus,  we  find  that  for  Run  1  u =  373.  In  fact, 
the  choice  of  s,n}  was  such  that  in  the  steady-state,  the  total  number  of 


particles  in  the  system  would  remain  close  to  the  initial  value,  a  choice 
which  was  not  fortuitous,  but  based  on  numerical  experimentation  with 
shorter  computer  runs. 

4.1  Evolution  of  the  Total  Number  of  Particles 

A  global  indication  of  the  behavior  of  Run  1,  as  it  evolved  through  the 
transient  phase  and  into  a  quasi-steady-state,  is  to  be  seen  in  the  time 
histories  of  the  total  electron  and  ion  populations,  Ne(t)  and  N,(t).  In 
Fig. (2),  we  plot  Ne(t)  and  Ni(t)  over  the  time  interval  0  <  a )ctt  <  1500.  In 
this  figure,  a  first  feature  to  be  noted  is  that  both  curves  show  a  sudden 
initial  drop.  This  sudden  loss  of  particles  results  from  the  scrape-off  of 
electrons  and  ions  over  a  layer  near  the  wall,  comparable  in  thickness  to 
the  gyro-radius  of  each  species.  The  electron  loss  (A Ne  &  —500)  occurs 
almost  instantaneously,  over  a  few  time  steps,  while  the  more  important  ion 
loss  (A Ni  «  —3500)  takes  place  over  the  longer  time-span  corresponding 
to  an  ion  rotation,  uc,t  «  27r.  The  ratio  of  the  numbers  of  particles  lost  is 
roughly  proportional  to  the  depth  of  the  scrape-offlayer,  with  AJV,/AiVe  ~ 
P,/fie  =  6.32. 

Immediately  after  this  initial  loss  due  to  scrape-off,  both  particle  popu¬ 
lations  begin  to  increase  in  response  to  the  external  source.  Note  that  in  the 
absence  of  any  outward  transport,  the  total  number  of  particles  of  either 
species  would  increase  roughly  five-fold  during  the  time-span  of  the  simu¬ 
lation,  since  the  filling-up  time  for  adding  15000  particles  is  Tjm  =  373u^“' . 
However,  this  is  not  the  case.  First,  the  initial  rates  of  increase  of  the  pop¬ 
ulations  are  lower  than  the  rate  of  injection  from  the  source.  Thus,  from 
Fig.(2)  we  find  that  dNe(t  rs  0 )/dt  rs  dNi(t  «  0 )/dt  =  0.67,  smaller  than 
=  1.005.  Furthermore,  both  dNe/dt  and  dNi/dt  decrease  with  time, 
with  a  noticeable  “bend”  occuring  in  the  population  curves  at  u)cit  k,  100. 
After  a  fairly  long  (0  <  ua^t  ~  1000),  but  not  very  pronounced  transient,  the 
system  settles  into  a  fluctuating  steady-state,  with  the  time-averaged  values 
of  dNe/dt  and  dNi/dt  near  zero.  In  the  steady-state,  we  have  Ne  «  16500 
and  Ni  ~  14000,  for  which  the  spatially-averaged  densities  are  such  that 
uFp,/wc,  Rs  1.11  and  5Tpe/u>ce  ~  0.19.  The  existence  of  the  near  steady-state 
in  the  presence  of  a  constant  source  of  particles  is  a  clear  indication  of 
the  existence  of  an  equal  and  compensating  outward  transport  of  particles, 


Figure  2:  Evolution  of  the  total  electron  and  ion  populations  in  Run  1;  the 
initial  values  are  Ne(0)  =  Ni( 0)  =  15000.  The  initial  scrape-off  of  particles, 
where  Ne  drops  to  14500  and  N,  to  12500.  occurs  for  ucit  <  2n  and  appears 
instantaneous  on  the  time-scale  of  the  figure. 

from  the  plasma  to  the  wall. 

4.2  Evolution  of  the  Electrostatic  Potential 

In  this  subsection,  we  shall  discuss  the  evolution  of  the  electrostatic  poten¬ 
tial  over  the  time  interval  0  <  u>clt  ~  1000.  This  interval  roughly  corre¬ 
sponds  to  the  the  transient  phase  of  Run  1  as  seen  in  Fig. (2). 

In  Figs. (3),  we  have  displayed  a  series  of  snapshots  of  <j>(x,  y,  t )  taken  at 
successive  moments  in  time.  Fig. (3a),  taken  at  ujcit  =  15,  shows  the  essen¬ 
tially  y-uniform  sheath  which  has  formed  at  the  beginning  of  the  evolution 
of  the  system,  after  21/2  rotations  of  the  ions.  This  sheath  is  due  to  the 
initial  loss  of  ions  which  have  impacted  into  the  wall,  resulting  in  a  layer  of 
depletion  of  positive  charge  over  a  depth  of  about  2  x  p;  ~  12,  and  a  net 
positive  charge  on  the  wall.  The  result  is  a  large  potential  drop  from  the 
wall  into  the  plasma,  which  in  Fig. (3a)  is  approximately  eA <t>/T,  %  —1.1. 
The  occurrence  of  a  potential  drop  from  wall  to  plasma  is  in  sharp  contrast 
to  the  situation  in  the  unmagnetized  sheath.  In  the  latter,  the  equilibrium 


configuration  is  dominated  by  the  electron  flow  to  the  walls,  and  this  results 
in  a  potential  rise  from  the  wall  into  the  plasma. 

The  sheath  shown  in  Fig. (3a)  contains  a  non-uniform  electric  field  Ex(x) 
pointing  into  the  plasma,  with  maximum  intensity  at  the  wall.  Because 
of  the  presence  of  the  external  magnetic  field,  this  electric  field  induces 
a  downward  E  x  B  drift  of  electrons  and  ions  parallel  to  the  wall,  with  a 
maximum  drift  velocity  which  is  very  close  to  the  ion  thermal  speed,  vy(x  = 
0)  =  —Ex(x  =  0)/Bo  «  —  vtt.  Because  the  electric  field  is  nonuniform,  the 
resulting  flow  is  strongly  sheared  in  velocity,  with  a  shear  length  of  order 
2  x  pi.  As  sheared  flow  is  in  general  unstable,  it  might  be  expected  that  the 
initial  structure  shown  in  Fig. (3a)  will  not  persist,  with  the  flow  vulnerable 
to  the  Kelvin- Helmholtz  instability  This  is  precisely  what  is  observed  in 
the  subsequent  evolution  of  the  potential:  y-dependent  ripples  (which  were 
already  visible  in  Fig. (3a)  )  grow  into  larger  perturbations,  which  are  both 
amplified  and  convected  by  the  E  x  B  flow.  This  process  appears  to  occur 
over  a  time-scale  (e-folding  time)  which  is  of  order  t  ~  20u»~1 .  In  Fig. (3b), 
at  Udt  =  62.5,  we  see  a  state  with  two  large  and  two  or  three  smaller 
perturbations  of  the  electrostatic  potential,  with  a  wavelength  which  we 
estimate  to  be  Ay  «  60.  A  measure  of  the  amplitudes  of  the  perturbations 
is  the  potential  drop  in  the  cross-sections  of  <f>(x,y)  taken  at  x  —  Lxj 2, 
which  range  from  e8(j>/Ti  ~  —0.08  to  e8<j>/T{  ~  —0.2. 

As  the  amplitudes  of  the  potential  perturbations  seen  in  Fig. (3b)  grow, 
they  undergo  a  first  saturation  by  forming  small  vortices  with  the  same 
characteristic  length,  Ay  ss  60.  These  vortices,  generated  in  the  time-span 
60  <  u>cit  ~  200,  are  short-lived  because  they  rapidly  coalesce  with  each 
other,  forming  a  smaller  number  of  vortices  with  Ay  ta  120.  There  is  also 
simultaneous  competition  from  the  mode  with  Ay  ss  120,  which,  as  we  shall 
see  in  the  next  section,  has  a  linear  growth  rate  comparable  to  that  of  the 
Ay  ~  60  mode.  The  overall  effect  is  to  favor  the  longer  wavelength.  As 
a  result,  Fig. (3c),  taken  at  wCJt  =  125,  shows  a  mixed  state,  :n  which  two 
(coalescing)  small  vortices  (Ay  ~  60)  coexist  with  a  large  vortex  (Ay  ss  120). 

Thus,  the  dynamics  of  the  system  favor  long  wavelengths  for  further 
amplification.  By  u ;c,t  =  250,  Fig. (3d),  only  two  large  vortices,  each  with 
e8d>/Tt  Rs  —1  and  Ay  ~  120,  have  survived.  These  vortices  drift  parallel 
to  the  wall,  in  accordance  with  the  local  E  x  B  drift,  and  can  be  consid¬ 
ered  quasi-stable,  because  they  maintain  a  constant  amplitude  and  survive 
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Figure  3:  Snapshots  of  the  electrostatic  potential  in  Run  1,  0  <  u jc,t  < 
437.5.  Left:  contour  plots  of  <j>(x,y);  Right  sections  of  <j>(x,y)  for  con¬ 
stant  x  =  16  =  Lx/2.  Note  that  in  the  contour  plots,  the  vertical  scale 
is  highly  compressed:  without  the  compression,  the  vortices  will  appear 
approximately  circular  ( continued  on  the  next  page). 
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Figure  3:  ( continued  from  the  previous  page)  Snapshots  of  the  electrostatic 
potential  in  Run  1,  0  <  uic,t  <  437.5.  Left :  contour  plots  of  <j>(x,y):  Right 
sections  of  <j>(x,y)  for  constant  x  —  16  =  Lx/2.  Note  that  in  the  contour 
plots,  the  vertical  scale  is  highly  compressed:  without  the  compression,  the 
vortices  will  appear  approximately  circular. 
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through  several  transits  across  the  system 

Let  us  clarify  the  connection  between  the  potential  structures  which 
we  have  called  “vortices’,  and  the  particle  motion.  If  we  neglect  finite 
gyroradii,  and  assume  that  the  particles  move  only  according  to  the  E  x  B 
drift  of  their  guiding  centers,  with  v  =  ExB  /  B2  =  z  x  V(j)/B ,  then  the 
electrostatic  potential  is  strictly  equivalent  to  a  two-dimensional  stream 
function  for  the  particle  flows  (this  assumption  is  better  for  the  electrons, 
because  of  their  small  gyroradius,  and  approximate  for  the  ions,  because 
their  gyroradius  is  a  sizeable  fraction  of  the  sheath  thickness).  Thus,  the 
closed  contours  of  the  potential  wells  tend  to  correspond  to  flows  with 
nonzero  circulation,  and  it  is  appropriate  to  refer  to  these  structures  as 
“vortices”. 

The  two- vortex  state  of  Fig. (3c)  is  in  fact  only  quasi-stable,  and  does 
not  last  beyond  u)c,t  «  400,  because  the  vortices  are  vulnerable  to  the  same 
coalescence  instability  which  engulfed  their  smaller  precursors.  Figs.(3d.e,f ) 
(ujc,t=  250,  375  and  437.5)  show  the  progression  from  the  state  with  two 
medium-amplitude  vortices  (e8<i>/T,  22  —1)  to  a  state  with  a  single,  large- 
amplitude  vortex  (eS'j/T,  «  —2.2).  The  process  appears  to  occur  with  one 
vortex  growing  larger  than  the  other,  with  the  vortices  then  attracting  each 
other  (Fig.(3e)),  and  finally,  with  the  smaller  vortex  merging  into  the  larger 
one  (Fig.(3f)). 

Fig. (4),  provides  a  more  detailed  view  of  vortex  formation  and  coa¬ 
lescence  than  is  available  from  the  individual  snapshots  of  Figs.(3).  We 
have  produced  a  time-series  of  cross-sections  of  the  electrostatic  potential 
<£(x,y,f),  taken  along  the  midplane  of  the  simulation  region,  at  a  fixed 
x  =  Lx/2  =  16.  These  cross-sections  are  identical  to  those  shown  in  the 
right-hand  column  of  Figs. (3)  and  were  produced  by  sampling  the  poten¬ 
tial  every  100  time  steps,  up  to  u;cxt  =  1045.  The  resulting  information  is 
the  two-dimensional  function  <psec{y^i)  ( x  being  fixed  and  ignorable)  and 
in  Fig. (4),  we  have  displayed  a  contour  plot  of  this  function.  In  this  form, 
the  vortices  of  Figs. (3)  immediately  appear  as  valleys  in  the  contour  map 
of  4>,'C(y,t).  Note  that  as  the  system  is  periodic  in  y.  perturbations  exiting 
at  y  =  0  immediately  reappear  at  y  =  256.  The  structure  of  oblique  bands 
in  the  plot  make  clear  the  steady  negative  drift  velocity  of  the  vortices. 

In  Fig. (4),  the  potential  profile  4>[y)  can  be  determined  at  any  one  time 
by  drawing  a  horizontal  fine  through  the  plot.  For  instance,  a  horizontal 


line  drawn  at  ucii  =  250  intersects  the  tracks  of  two  vortices,  in  agreement 
with  the  picture  of  Fig. (3d).  If  we  now  sweep  this  horizontal  line  upward, 
we  can  follow  in  detail  the  coalescence  process  which  was  only  outlined  in 
Figs.(3d-f).  Near  u>c,t  =  330,  the  left-hand  and  larger  vortex  undergoes  a 
slight  deceleration  (its  y-directed  velocity  becomes  less  negative);  somewhat 
later,  at  uc,t  =  355,  the  right-hand  and  smaller  vortex  undergoes  a  more 
sizeable  acceleration  (its  y-directed  velocity  becomes  more  negative).  There 
follows  a  period  of  “coasting”,  375  <  <  425,  where  the  two  vortices  are 

in  close  proximity.  Finally,  by  u >cit  =  425,  the  two  vortices  merge  together,  a 
process  in  which,  because  of  the  disparity  in  sizes,  the  larger  vortex  appears 
to  engulf  the  smaller  one. 

Fig. (4)  provides  quantitative  information  on  the  drift  velocities  of  the 
vortices,  the  instantaneous  drift  velocity  being  equal  to  the  reciprocal  of 
the  slope  of  the  vortex  trajectory  in  ( y,t )  space.  A  notable  feature  of  these 
velocities  is  that  the  average  vortex  speed  (averaged  over  one  transit  of  the 
system)  is  almost  constant.  For  instance,  at  u c,t  =  150,  when  the  vortex 
amplitudes  are  of  order  e8<p/Ti  as  —0.5,  we  find  a  vortex  drfit  velocity 
vy  =  —0.069  (\vy\/vti  —  0.44).  At  u >c,t  =  1000,  when  the  vortex  amplitude 
is  e84>/T,  ss  —2,  we  find  a  comparable  drift  velocity,  vy  —  —0.075  (|ey|  /vtl  = 
0.47). 

We  also  note  in  connection  with  Fig. (4)  that  in  the  final  quasi-steady 
state,  u:cit  >  400,  the  shear  layer  is  still  locally  unstable  and  is  capable  of 
amplifying  perturbations  at  some  distance  from  the  main  vortex.  Thus, 
smaller  “satellite”  vortices  are  produced.  For  instance,  in  Fig.  ( 4 )  we  see 
the  appearance  at  950  of  a  small  vortex  at  a  distance  ly  s;  90  to 

the  left  of  the  main  vortex.  At  u>Cit  ss  1020,  the  satellite  vortex  reaches 
a  maximum  depth  of  eS<t>/T,  «  -0.6,  or  about  1/3  the  depth  of  the  main 
vortex.  It  is  eventually  absorbed  by  the  main  vortex,  by  w 'Clt  =  1050. 

Finally,  we  note  that  as  the  vortices  grow  and  coalesce,  the  shear  layer 
undergoes  a  concommitant  broadening.  In  Figs.(5a.b),  we  show  the  pro¬ 
files  for  the  y-averaged  electric  field  Ex(x),  and  for  the  y-averaged  charge 
density,  p(x)  =  dEx(x)/dx,  at  three  successive  times.  u-'c>t  =  15.  250,  1075. 
Since  t’y  =  — Ex/ B ,  Fig. (5a)  is  also  a  plot  of  the  profile  of  velocity  in  the 
shear  layer.  Fig. (5b)  shows  that  the  inflection  point  x0  of  the  velocity 
profile  ( i'y(xo)  =  —  jjdp(x0)/ dx  =  0)  moves  outward,  and  that  the  overall 
profiles  broaden.  Note  that  at  all  times  the  average  shear  profile  is  linearly 
unstable,  with  the  inflection  point  and  regions  of  opposite  curvature  needed 
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Figure  5:  Cross-sections  of  the  sheath  in  Run  1.  for  u ictt  =  15,  250  and  1075; 
(a)  {/-averaged  electric  field  Ez{x)\  (b)  y-averaged  charge  density  p{x). 

for  the  Kelvin-Helmholtz  instabihty.  This  feature  shows  that  the  satura¬ 
tion  mechanism  is  not  quasilinear  (relaxation  of  the  space-averaged  shear 
profile),  but  indeed  strongly  nonlinear. 

4.3  Evolution  of  the  Fourier  Modes 

We  now  complement  the  observations  of  the  previous  paragraphs,  in  which 
we  identified  potential  structures  in  real  {x.y)  space,  with  a  discussion  of  the 
evolution  of  the  Fourier  transform  of  the  potential.  o{x,ky,t).  In  Figs.  (6), 
we  show  the  temporal  behavior  of  the  first  four  Fourier  modes  with  ky  ^  0. 
with  the  Fourier  transforms  taken  along  the  midplane  of  the  system,  at  a 
fixed  x  —  Lzj 2  =  16.  We  have  plotted  both  real  and  imaginary  parts  of 
o(j",  ky ,  t).  The  time  interval  is  identical  to  that  of  Fig. (4),  0  <  u,-cit  <  1050. 

A  first  salient  feature  of  Figs. (6)  is  that  the  phase  of  exponential  growth 
of  any  one  mode  is  quite  limited  in  time.  A  magnification  of  the  beginning 
of  Figs. (6)  shows  that,  starting  from  noise,  the  exponential  growth  of  the 
modes  does  not  extend  beyond  roughly  two  e-foldings.  Thus,  the  progres¬ 
sion  of  the  system  toward  the  steady-state  occurs  for  the  most  part  with 
what  appears  to  be  algebraic  growth.  This  behavior  is  consistent  with  the 


Figure  6:  Temporal  dependence  of  the  Fourier  modes  in  Run  1.  modes 
m  =  1,  2,  3  and  4.  The  Fourier  transform  is  taken  along  the  midplane  of 
the  system,  at  a  fixed  x  =  Lx/2  =  16.  Solid  lines:  real  part;  dashed  lines: 
imaginary  part. 


description  of  the  evolution  of  <£(x,y,  f)  given  in  the  previous  section,  in 
which  we  saw  that  a  first  nonlinear  state,  with  formation  of  small  vortices, 
occured  as  early  as  ioc,t  ss  100  (Fig. (3c).  Beyond  this  point  in  time,  the 
dominant  mechanisms  for  the  evolution  of  the  instability  are  nonlinear, 
consisting  of  the  continued,  nonlinear  growth  of  individual  vortices,  and 
of  their  pairwise  coalescence.  A  qualitative  analogy  which  comes  to  mind, 
is  the  nonlinear  evolution  of  MHD  tearing- modes,  in  which  the  algebraic 
growth  and  the  coalescence  of  magnetic  islands  bear  some  resemblance  to 
the  behavior  of  our  electrostatic  vortices[l9]. 

The  coalescence  of  the  smaller  vortices  into  larger  ones  is  a  process 
indicative,  in  wavenumber  space,  of  an  inverse  cascade.  In  Figs. (6a. b).  we 
have  a  clear  indication  of  the  twn-to-one  vortex  coalescence  which  was  seen 
to  occur  in  real  space  in  the  time-interval  250  <  u jcxi  <  437  (Figs.(3d-f )). 
the  two- vortex  state  corresponding  to  the  m  =  2  mode,  and  the  one- vortex 
state  to  the  m  =  1  mode.  The  process  begins  with  the  sudden  excitation 
of  the  fundamental  mode  m  =  1,  at  about  ujc,t  =  330,  coinciding  with 
the  beginning  of  spatial  convergence  of  the  two  vortices  which  was  seen  in 
Fig. (4)  The  actual  coalescence  process,  m  which  one  vortex  engulfs  the 
other,  occurs  at  u ic,t  ss  400,  and  is  signaled  by  a  sizeable  perturbation  of 
the  m  =  2  mode. 

In  Fig. (7)  we  show  the  relative  amplitudes  of  the  Fourier  modes  at  u jcit  = 
1050,  at  which  time  the  modes  have  essentially  settled  into  a  quasi-steady- 
state.  This  figure  shows  that  the  first  5  or  6  modes  are  dominant,  with 
the  overall  Fourier  spectrum  strongly  cutoff  beyond  kvp,  —  1  Finally,  let 
us  note  that  the  short- wavelength  modes  not  only  have  smaller  amplitudes, 
but  axe  also  more  strongly  fluctuating,  as  can  be  seen  by  a  casual  inspection 
of  Fig. (6)  This  visual  impression  is  confirmed  in  the  analysis  of  the  steady- 
state  power  spectra  P(*\  ky),  which  we  present  in  Section  6.5. 

4.4  Minimum  Length  for  a  Two- Vortex  Steady-State 

In  the  present  simulation,  we  did  not  determine  the  minimum  length  in  y 
required  for  the  existence,  in  the  steady-state,  of  more  than  one  large  vortex. 
Furthermore,  limitations  on  computational  ressources  did  not  permit  us  to 
explore  longer  systems  in  which  such  a  state  might  appear.  We  would  like 
to  argue  however,  using  qualitative  arguments,  that  a  minimum  length  for  a 
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Figure  7:  Relative  amplitudes  of  the  Fourier  modes  |d(ktf,  x.  <)|  at  t  = 
1050.  Run  1.  for  x  =  Lxj 2  =  16. 

two-vortex  steady-state  must  exist.  In  other  words,  we  are  suggesting  that 
it  is  physically  reasonable  to  suppose  that  vortex  coalescence  can  go  not 
on  indefinitely,  and  that  the  final  state  of  any  sufficiently  long  system  will 
consist  of  more  than  one  large-amplitude  vortex.  To  support  our  assertion, 
let  us  consider  a  scenario  where  such  a  single- vortex  state  might  initially 
prevail.  We  then  argue  that  it  may  nonetheless  result  in  a  many-vortex 
state,  provided  that  the  system  is  long  enough. 

Let  us  assume  that  the  single,  final  vortex  is  roughly  circular,  and  oc¬ 
cupies  at  any  one  time  a  limited  extent  in  y,  comparable  to  the  sheath 
thickness  in  x.  that  is  a  few  ion  gyroradii.  Then,  we  cannot  overrule  the 
possibility  that  in  a  very  long  system  the  shear  layer  may  produce  satel¬ 
lite  vortices  at  some  distance  from  the  main  vortex,  as  was  seen  in  Fig. (4) 
at  u :c,t  %  1000.  Provided  that  the  system  is  long  enough,  and  that  the 
secondary  vortices  are  appearing  far  enough  from  the  main  vortex,  these 
satellite  vortices  should  have  the  time  to  develop  large  amplitudes  of  their 
own.  even  if  they  eventually  merge  with  other  vortices  in  the  system.  Thus, 
in  the  steady-state  the  system  may  contain  many  vortices.  In  connection 
with  this  scenario,  we  would  also  like  to  argue  that  the  oldest  and  largest 
vortices  may  not  grow  indefinitely  by  engulfing  the  satellite  vortices.  Pre- 


sumably,  very  large  amplitude  vortices  with  \ef>o/T,  I  »  1  are  unstable,  and 
can  shed  energy  either  by  interaction  with  the  wall,  or  perhaps  emission  of 
waves.  The  fined  result  will  be  a  self-limiting  state  consisting  of  many  large 
(but  not  indefinitely  large)  amplitude  vortices,  in  constant  interaction  with 
each  other. 

The  qualitative  arguments  outlined  above  are  clearly  insufficient,  and 
should  be  supported  by  future  numerical  simulation  of  longer  systems. 

4-5  Saturation  Mechanisms:  Climax  and  Coalescence 

We  shall  now  make  a  comparison  between  the  transient  plasma  behavior 
observed  in  Run  1.  and  the  transient  behavior  of  the  fluid  Kelvin- Helmholtz 
instability,  as  observed  m  purely  hydrodynamic  fluid  simulations  Our  aim 
is  to  provide,  through  this  analogy,  qualitative  explanations  of  the  satura¬ 
tion  and  coalescence  phenomena  seen  in  Run  1.  The  basis  of  the  analogy 
is  the  close  similarity,  in  the  fluid  approximation,  of  the  plasma  cross -field 
equations,  to  be  derived  in  Section  5,  with  the  hydrodynamic,  Navier-Stokes 
equation  (21). 

We  note  here  that  a  rigorous,  analytic  solution  of  the  Navier-Stokes 
equations  for  the  nonlinear  evolution  of  the  fluid  Kelvin-Helmholtz  insta 
bility  is  so  far  lacking.  Though  a  variety  of  special  steady-state  solutions  can 
be  found  (periodic,  soliton-like,  localized),  there  is  no  rigorous  procedure  to 
predict  the  appearance  of  a  given  kind  of  nonlinear  structure  from  arbitrary 
initial  conditions  (as  can  be  done,  say,  in  the  case  of  the  Ivorteweg-deVries 
equation).  Because  of  this  state  of  affairs,  we  expect  that  a  description 
of  the  nonlinear  dynamics  of  our  system,  even  in  the  fluid  approximation, 
must  remain  at  best  semi-quantitative. 

As  we  saw,  the  general  trend  in  Run  1  was  the  following:  after  a  few 
e-folding  times  the  linear  instability  saturated  into  a  state  consisting  of  sev¬ 
eral  quasi-stable  vortices,  which  eventually  coalesced  into  larger  structures, 
leading  to  a  final  unique  and  stable  vortex.  This  evolution  is  quite  similar 
to  that  observed  in  numerical  simulations  of  two-dimensional,  sheared  fluid 
flows.  For  the  purpose  of  the  present  analogy,  we  shall  especially  refer  to 
the  work  of  Corcos[20,21,22]. 

Corcos[20,21,22j  has  presented  a  qualitative  model  for  the  saturation  of 
single  vortices  in  the  fluid  Kelvin-Helmholtz  instability,  starting  from  small 
amplitudes.  The  saturated  state  has  been  termed  the  climax  state.  We 


believe  that  this  model  is  qualitatively  valid  for  the  plasma  instability  as 
well,  and  we  shall  sketch  its  outlines.  Let  us  first  emphasize  the  connection 
between  vorticity  and  charge  density  in  the  plasma  fluid.  If  the  plasma 
flows  are  given  by  v  =  E  x  B/fl,  then  the  vorticity  is  simply: 

uz  =  (V  x  v),  =  -p/Be0  (1) 

where  p  is  the  local  charge  density.  Thus,  at  any  time  we  can  identify 
vorticity  with  (minus)  the  charge  density. 

The  Corcos  model  for  the  evolution  to  the  “climax”  state  depends  on  the 
reciprocal  interaction  of  the  vorticity  and  the  fluid  flow,  with  the  vorticity 
driving  the  flow  through  Ampere’s  law,  and  being  itself  advected  by  the 
flow  through  the  continuity  equation.  The  model  proceeds  as  follows:  as 
the  initially  small-amplitude  fields  grow,  they  divide  the  fluid  into  free  and 
trapped  flows,  the  trapped  flows  corresponding  to  the  growing  vortices. 
Now,  the  Kelvin-Helmholtz  instability  grows  by  extracting  vorticity  from 
the  free-flowing  shear  layer.  This  vorticity  is  convected  into  the  trapped 
regions,  where,  once  trapped,  it  becomes  unavailable  as  a  source  of  free- 
energy.  The  convection  of  vorticity  from  free  to  trapped  regions  proceeds 
in  two  steps;  first,  the  vorticity,  which  is  initially  diffuse  throughout  the 
shear  layer,  is  concentrated  by  the  large-amplitude  flows  into  a  narrow  band 
along  the  vortex  separatrix;  this  band  of  vorticity  is  then  advected  into  the 
vortex  core  at  a  point  roughly  at  the  end  of  the  separatrix  (at  the  point  of 
maximum  excursion  into  the  shear  layer),  by  a  combination  of  stretching 
and  rolling  of  the  stream-line  defining  the  separatrix.  Because  the  total 
supply  of  vorticity  is  constant  in  this  inviscid  flow,  the  instability  will  cease 
to  grow  when  all  of  the  vorticity  has  been  extracted  from  the  shear  layer 
and  concentrated  inside  the  vortices  (because  of  the  uncertainties  of  the 
model,  in  fact  it  only  requires  that  an  an  appreciable  fraction  of  the  total 
vorticity  be  trapped). 

Let  Qvot  denote  the  net  charge  trapped  inside  the  vortex,  and  Q0  the 
initial,  net  charge  in  one  wavelength  of  the  Kelvin-Helmholtz  instability. 
Then,  on  account  of  Eq.(l),  the  climax  condition  for  the  saturation  of  the 
individual  vortices  can  also  be  expressed  as: 


Qvor  «  Qo  (2) 

In  other  words,  the  vortex  saturates  when  it  has  trapped  most  of  the  net 


charge  resident  in  the  shear  layer.  This  occurs  when  the  separatrix  width 
becomes  comparable  to  the  width  of  the  shear  layer,  and  is  indeed  what 
is  observed  in  Figs. (3).  Another,  semi-quantitative  consequence  of  Corcos’ 
model  is  that  the  growth  of  the  vortex  should  proceed  at  an  algebraic  rate 
throughout  the  finite-amplitude  part  of  its  evolution,  a  feature  which  was 
observed  in  the  graphs  of  Figs. (6). 

Let  us  now  comment  on  the  process  of  vortex  coalescence,  a  process 
which  is  also  a  generic  feature  of  the  hydrodynamic  simulations,  and  which 
immediately  follows  the  climax  of  individual  vortices.  The  vortices  coalesce 
as  a  result  of  a  nonlinear  instability,  the  so-called  pairing  instability  first 
derived  by  Lamb  for  point  vortices[23]  This  instability  has  since  then  been 
analysed  for  more  realistic  vortex  configurations[24]. 

For  a  row  of  two-dimensional  point-vortices,  initially  equally  spaced  at  a 
distance  6,  and  each  with  net  circulation  ruor  =  /  u>*  dxdy ,  one  finds  that  the 
arrangement  is  unstable,  with  a  maximum  Unear  growth  rate  corresponding 
to  a  motion  in  which  the  vortices  converge  in  pairs,  the  pairing  instability. 
The  growth  rate  for  this  mode  is[23]i 

7pa.r  =  “T  (3) 

In  Section  5.2,  we  shall  see  that  for  a  sheax  layer  of  width  6,  with  charac¬ 
teristic  velocity  the  maximum  Unear  growth  rate  scales  as: 

1L  ~  ~r  f°r  b  (4) 

b 

where  we  have  ignored  all  numerical  factors.  To  compare  Eqs.(3)  and  (4), 
we  note  that  the  total  circulation  To  in  one  wavelength  of  the  initial  shear 
layer  is: 

To  =  bV0  (5) 

Using  Eqs.(3),  (4)  and  (5),  we  find  that: 


Now  in  Corcos'  model,  the  vortices  saturate  precisely  when  rior/r0  ~  1. 
This  imphes  that  as  soon  as  the  vortices  saturate,  the  pairing  instability  will 
immediately  set  in,  with  a  growth  rate  of  order  ->lir  —  -  /  In  agreement 


with  Figs. (3),  this  precludes  a  long-lived  many-vortex  state,  at  least  in  the 
early  phases  of  evolution,  when  the  vortices  are  closely  packed  in  a  periodic 
or  near-periodic  fashion. 

4.6  Behavior  when  k\\  0 

We  briefly  explored  the  behavior  of  the  sheath  when  k\\  ±  0.  In  our  two- 
dimensional  simulations,  this  was  done  by  tilting  the  magnetic  field  at  a 
finite  angle  to  the  z-axis.  We  considered  a  system  half  as  large  as  in  Run  1, 
with  Ly  =  128  and  with  half  the  particle  injection  rate,  but  with  otherwise 
identical  parameters. 

In  a  first  configuration,  we  tilted  the  magnetic  field  in  the  direction 
of  the  y-axis,  by  setting  B  =  (0,  B  sin  a,  B  cos  a),  with  an  angle  a  = 
j(me/m,)1/2  =  4.53°.  In  this  configuration  the  magnetic  field  remained 
parallel  to  the  wall,  but  allowed  particle  flow  along  the  field  fines  in  the 
y  direction,  with  a  parallel  velocity  of  order  t>||  ~  utsina,  where  vt  is  the 
thermal  velocity. 

We  expected  important  changes  in  the  sheath  dynamics  to  occur  when 
the  electrons  could  flow  at  a  parallel  velocity  comparable  to  their  E  x  B 
drift  velocity.  Since,  as  we  shall  see,  the  E  x  B  velocities  in  the  sheath 
are  of  order  vp  ~  vtj,  we  expected  the  transition  to  occur  when  i>e||  ~  vn, 
or  at  a  “critical”  angle  a  ~  (me/m,j^2.  In  fact,  even  the  smaller  an¬ 
gle  a  =  i(me/m,)1/2  was  sufficient  to  greatly  modify  the  course  of  the 
simulation,  by  completely  suppressing  the  formation  of  vortices  and  edge 
turbulence.  While  the  system  generated  a  one- dimensional  sheath  almost 
identical  to  the  one-dimensional  sheath  generated  in  the  early  phases  of 
Run  1  (with  potential  drop  eA(f>/Ti  =  —1.4),  this  one-dimensional  sheath 
remained  absolutely  stable  in  time,  and  did  not  induce  particle  transport 
(most  particles  injected  into  the  system  remained  confined).  We  suspect 
that  this  stability  is  a  consequence  of  the  parallel  electron  motion,  which 
stabilizes  the  Kelvin-Helmholtz  instability  by  short-circuiting  potential  per¬ 
turbations  which  might  arise  in  the  shear  layer. 

In  a  second  configuration,  we  tilted  the  magnetic  field  in  the  direction  of 
the  x  axis,  thereby  making  the  field  lines  impinge  onto  the  walls.  This  was 
done  by  setting  B  =  (B  sin  /?,0,  B  cos  /?),  with  13  =  (me/m,Y^2  / 8  =  1.13°. 
Despite  the  small  value  of  the  tilting  angle,  this  configuration  resulted  in  a 
completely  different  type  of  sheath,  in  which  the  electron  flow  dominated 


the  sheath  dynamics.  The  wall  took  on  a  negative  surface  charge,  and 
the  electrostatic  potential  now  exhibited  a  rise  going  from  the  wall  to  the 
plasma,  of  order  eA<f>/T  ~  0.3.  Thus,  in  this  configuration,  the  cross-field 
sheath  was  suddenly  made  to  ressemble  the  unmagnetized  sheath,  and  we 
can  regard  the  introduction  of  the  tilting  angle  j3  as  a  “singular”  pertur¬ 
bation.  We  did  not  perform  a  systematic  study  of  the  sheath  with  ^  0; 
however  we  expect  that  it  is  well  described  by  the  work  of  Chodura[8]. 

The  results  of  the  two  cases  outlined  above  emphasize  that  the  sheath 
model  that  we  are  studying  in  this  paper  applies  strictly  to  flute-like  modes, 
with  a  very  small  feu.  For  our  k\\  —  0  model  to  be  appoximately  valid  when 
&II  7^  0,  requires  that  in  the  characteristic  time  needed  for  the  Kelvin- 
Helmholtz  instability  to  develop,  the  electrons  sample  only  a  fraction  of  a 
parallel  wavelength  Ay  =  2n/k\\.  This  imposes  the  condition  A||/ufe  »  7k//' 
where  'ykh  is  the  growth-rate  of  the  Kelvin-Helmholtz  instability.  As  we 
shall  see  in  Section  5.2,  7^  >  20a;"1  (Fig. (12)).  The  result  is  the  condition: 

A»>20t)  “•  (7) 

or  in  terms  of  parallel  wavenumber: 
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5  Fluid  Theory 

In  this  section,  we  derive  the  fluid  equations  for  the  cross-field  dynamics  of 
the  electrons  and  ions.  Our  aim  is  twofold:  in  the  present  section  we  use  the 
cross-field  equations  in  linearized  form,  and  predict  the  initial  growth  rates 
for  the  Kelvin-Helmholtz  instability.  In  Section  6  we  shall  then  consider 
nonlinear  solutions  of  these  equations,  and  attempt  to  fit  these  solutions  to 
the  final  steady-state  observed  in  Run  1.  Our  approach  is  similar  to  that 
taken  by  several  authors[l6,25). 

The  cross-field  equations  are  two  coupled,  nonlinear  evolution  equa¬ 
tions,  linking  the  electrostatic  potential  to  the  particle  density.  As  we 
shall  see,  there  is  a  close  analogy  between  these  plasma  equations  and  the 
two-dimensional,  inviscid  Navier-Stokes  equation[26].  This  analogy  was  the 
basis  for  the  discussion  of  Section  4.5,  in  which  we  compared  the  evolution 
of  the  plasma  simulation  to  that  of  hydrodynamic  simulations. 

Following  a  general  derivation  of  the  fluid  equations,  we  derive  the  lin¬ 
earized  form  of  the  cross-field  equations,  assuming  that  the  equilibrium 
conditions  at  the  edge  are  known.  To  determine  these  equilibrium  condi¬ 
tions,  we  fit  analytic  (tanh(i))  profiles  to  the  initial  velocity  and  density 
profiles  measured  in  Run  1.  We  then  solve  numerically  the  linearized  fluid 
equations,  and  obtain  estimates  for  the  growth  rates  and  frequencies,  and 
the  eigenfunctions  of  the  unstable  Kelvin-Helmholtz  modes.  This  Unear 
analysis,  which  assumes  a  nonuniform  density  profile  and  the  presence  of 
a  finite  conducting  boundary,  is  considerably  more  general  than  the  usual 
analysis  of  the  Kelvin-Helmholtz  instability[12],  which  is  done  for  a  constant 
density  profile  and  boundary  conditions  at  infinity.  The  results  of  our  linear 
analysis  are  then  compared  to  the  numerical  results  of  the  two-dimensional 
simulations,  with  the  aim  of  confirming  that  the  instabilities  seen  in  the 
simulations  can  indeed  by  ascribed  to  the  Kelvin-Helmholtz  instability  of 
the  edge  shear  layer.  We  find  reasonable  agreement  for  the  long-wavelength 
modes,  but  a  sharper  cutoff  for  the  shorter  wavelength  modes,  indicating 
strong  finite-larmor-radius  effects  not  accounted  for  in  the  fluid  theory. 

The  fluid  theory  assumes  that  we  are  modelling  low-frequency,  long- 
wavelength  phenomena,  with  |u;[  u,’ct  and  [k|p,  1.  Fortunately,  inspec¬ 

tion  of  the  self-consistent  simulation  results  shows  that  these  conditions  are 
indeed  approximately  satisfied.  Fig. ( 7 )  showed  that  the  dominant  Fourier 
modes  were  those  with  | |  p,  <  1,  and  an  examination  of  the  power  spec- 
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tra  P(u>,ky)  (Fig. (20),  to  be  discussed  in  Section  6.5)  shows  that  for  the 
long- wavelength  modes,  power  is  concentrated  at  frequencies  |u>|  <C  u,’c>. 

As  an  aside,  we  note  that  the  dominance  of  low-frequency,  long- wavelength 
modes  over  short-wavelength,  high-frequency  fluctuations  can  be  readily 
seen  in  a  computer  animation  of  the  potential  surface  <p(x,y,t),  which  we 
generated  from  data  produced  in  a  simulation  with  parameters  similar  to 
those  of  Run  l[5]. 

Thus,  the  fluid  theory  correctly  describes  the  dominant,  long-wavelength, 
low-frequency  features  of  the  cross-field  dynamics  (the  linear  instability,  the 
formation  of  vortices,  vortex  coalescence  and  vortex  structure),  but  will  not 
account  for  the  small- amplitude,  fluctuating  behavior  of  the  fields  at  short 
wavelengths,  kyp;  ~  1.  Presumably,  the  short-wavelength  phenomena  are 
correctly  described  only  by  the  full  kinetic  equations,  the  solution  of  which 
is  a  topic  for  future  research. 


5.1  The  Nonlinear  Cross-Field  Equations 

Our  derivation  is  similar  to  the  one  presented  by  Horton  et  aJ .  1 16)  In  the 
regime  u;  u>cl,  we  simphfy  the  electron  and  ion  momentum  equations  by 
approximating  the  electron  motion  by  the  E  x  B  drift,  and  the  ion  motion 
by  the  E  x  B  and  polarization  drifts: 


ve  =  -(E  x  z)  =  vE 


1  1 
v,  =  —  (E  x  z)  + 


3 


dt 


E  =  vE  +  vP 


In  Eq.(10),  d/dt  is  the  total  derivative, 
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With  the  electron  and  ion  flows  given  by  Eqs.(O.lO),  the  densities  must 
satisfy  the  continuity  equations: 
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—  ne  4- V  ■  (v£nc)  =0 
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—  n,  +  V  ■  ((V£- +  vP)ne)  =0 


(12) 


(13) 


Subtracting  Eq.(12)  from  Eq.(13),  we  obtain  a  transport  equation  for  the 
net  charge  separation: 


— (n,  -  ne)  +  V  •  (v£(n,  -  ne)  +  vPn,)  =  0 
Using  V  •  V£  =  0,  Eq.(14)  can  be  written: 

s(','-n')+iv(',4E)=o 

If  we  now  use  Poisson's  equation: 


(14) 


(15) 


V2<£  =  -~{n,  -  ne)  (16) 

Co 

and  substitute  for  n,  —  ne  in  Eq.(15),  we  obtain  a  single  equation  tying 
the  electrostatic  potential  to  the  ion  density.  To  this  equation,  we  append 
the  electron  continuity  equation  (12).  Using  the  quasineutrality  condition. 
ne  «  n,  E  n,  we  obtain  the  coupled  nonlinear  equations: 


m,e pu;2  d 
e 2  dt 


V20  + V- 


„-v,  =° 


tm 


7<n  = 0  (18) 

This  set  of  equations  is  similar  to  the  cross-field  equations  of  Horton  et 
al. [27] ,  differing  by  the  addition  of  the  first  term  in  Eq.(17),  which  allows 
for  a  finite  value  of  uct / upi .  Eqs.(17,18)  have  some  simpler  limiting  forms. 
^  ^  ^c.’  then  we  can  neglect  the  first  term  in  Eq.(17)  and  obtain: 


v(4S=0  (19) 

Ttn  =  °  <20> 

If  we  further  assume  that  n  is  everywhere  constant,  so  that  the  transport 
equation  (20)  is  trivially  satisfied,  we  obtain  a  single  nonlinear  equation  for 
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analytic  expressions  of  the  form: 

E. 

r0(x)  =  Vq  f0(  x ) 

(22) 

nu(x)  =  .V0  Uq(x) 

(23) 
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^+i(ixv«  V)  w  =  0 


(21) 


Eq.(21)  is  identical  to  the  two-dimensional  Navier-Stokes  equation,  with 
stream- function  0.  We  note  that  to  derive  Eq.(21)  from  Eq. ( 17)  does  not 
require  the  assumption  of  u »  u but  only  that  of  constant  density.  Fi¬ 
nally,  if  <C  u>£,  Eq.(17)  once  again  imphes  Eq.(21),  this  time  irrespective 
of  any  assumptions  on  the  behavior  of  the  density  profile.  Thus  Eqs.(  17.18) 
are  close  in  nature  to  the  two-dimensional  Navier-Stokes  equation,  Eq.(21), 
in  a  manner  that  is  not  very  sensitive  to  the  value  of  u^/ui^,. 


5.2 


Growth-Rates  of  the  Kelvin-Helmholtz  Instabil¬ 
ity 


We  shall  now  estimate  the  growth  rates  predicted  by  Eqs.(l7,18)  in  the 
presence  of  specified  zero-order  velocity  and  density  profiles.  In  what  fol¬ 
lows.  we  did  not  use  the  exact  numerical  profiles  obtained  in  Run  1.  but 
rather,  we  fitted  approximate  analytic  forms  to  these  profiles:  this  proce¬ 
dure  leads  to  a  fully  analytic  treatment  of  the  instability,  and  avoids  the 
problems  associated  with  finite  noise  in  the  numerical  results. 

In  Fig. (8),  we  show  a  snapshot  of  the  “exact”,  zero-order  profiles  as  they 
were  obtained  from  the  self-consistent,  two-dimensional  calculation  of  Run 
1  The  snapshot  is  taken  at  an  early  tune  in  the  evolution  of  the  system. 
ojcxt  —  15,  at  which  nine  the  sheath  is  still  uniform  in  y  To  predict  the 
growth  rates,  we  model  these  numerical  velocity  and  density  profiles  by 


where  we  have  defined  c0y(x)  =  v0(x).  In  Eqs. (22.23)  1’0  and  N0  are  con¬ 
venient  reference  velocities  and  densities.  We  have  chosen  the  functional 
dependences: 


c0(  r )  =  tanh(r/a  -  1)  -  1 


(24) 
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With  Eqs. (22-25),  the  equilibrium  density  far  from  the  wall  is  given  by  iV0, 
the  drift  velocity  at  the  wall  is  given  by  —(1  4-  tanh(l))  =  — 1.761  Vq,  d  is 
the  scale  length  for  the  density  gradient,  and  finally,  a  locates  the  position 
of  the  inflection  point  in  the  velocity  profile  ( p'(a )  =  0).  Fitting  to  the 
simulation  curves  of  Fig. (8),  we  find  the  values  (in  the  units  of  ES2): 


Vq  =  0.0812  a  =  4.0  (26) 

d  =  8.0  u£/u£(n  =  N0)  =  0.73  (27) 

or.  in  physical  parameters: 

V0/vtt  =  0.514  a/ Pi  =  0.63  (28) 

d/p,  =  1.26  (29) 

The  analytic  profiles  which  result  from  this  choice  of  parameters  are  shown 
in  Fig. ( 9 ) .  It  is  convenient  do  introduce  the  parameter  a,  defined  as: 
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“ !p,in  =  iV°) 

Let  us  now  linearize  Eqs. (17, 18),  with  the  equilibrium  given  by  Eqs. (22. 23). 
Writing  O  =  <i>o(x)  +  <t>i(x,y,t)  and  n  =  n0(x)  +  rii(x,y, <),  we  find  after 
some  algebra: 

f  d  9  \  (  2  d'0(x)  d  ,  \ 

‘  (,,;:(x)  +  =  0  (31) 
In  this  equation,  we  have  not  yet  factored-out  VQ  from  vQ(x).  A  simplifi¬ 
cation  occurs  in  deriving  Eq.(31),  in  that  does  not  appear  in  Eq.(17). 
Thus,  to  first  order  Eqs. (17)  and  (18)  are  decoupled,  and  one  need  linearize 
only  Eq.(  17)  to  obtain  a  single  equation  for  Oi-  as  we  have  done  in  (31). 
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Figure  8:  Velocity  and  density  profiles  for  Run  1,  at  uclt  =  15;  (a)  y- 
averaged  velocity  vy(x)  =  v0(x )  and  charge-density  profiles  (full  and  dashed 
hnes  respectively);  note  that  the  minimum  of  p(x ),  at  x  =  a,  determines 
the  inflection  point  v "  =  0;  (b)  density  profiles  for  electrons  and  ions;  the 
zero-order  density  in  Eq.(23)  is  n0{x )  =  n,(x),  and  d  is  the  scale  length  in 
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Figure  9:  Shapes  of  the  analytic  velocity  and  density  profiles  used  in  pre¬ 
dicting  the  growth-rates  of  the  Kelvin- Helmholtz  instability. 

A  rather  lengthy  eigenvalue  equation  is  obtained  by  assuming  an  exp(i(kyy 
uit))  dependence  for  the  normal  modes: 


K(x)  d  ,  ,  f  K 


a  +  h0(x)dx  \w  —  kyv0(x) 


«o(xK(*) 


a  +  h0(x) 


Eq.(32)  can  be  made  to  appear  more  familiar,  if  we  realize  that  for  a  con¬ 
stant  density  profile,  h'0(x)  =  0,  or  in  the  limit  of  low  peak  densities,  a  — »  oo, 
Eq.(32)  reduces  to  the  well-known  Rayleigh  stability  equation[28,29].  The 
new  terms  in  Eq.(32)  are  those  proportional  to  the  density  gradient  h’0(x). 

We  now  consider  the  numerical  solution  of  Eq.(32),  for  arbitrary  density 
profiles.  To  find  the  frequencies  and  growth  rates  of  the  unstable  modes, 
it  is  convenient  to  solve  Eq.(32)  not  as  an  eigenvalue  problem,  but  directly 
in  the  time-dependent  form  (31),  as  an  initial-value  problem[15].  In  this 
approach,  we  introduce  an  intermediate  quantity  defined  as  the  term  in 
Eq.(31)  that  is  operated  upon  by  the  convective  derivative  ( dt  +  v0(x)Oy). 
Eq.(31)  is  then  re-written  as  a  coupled  system: 


a  +  hQ(x)  dx 


•'y  I  ^1  -  Vl 


Eqs.(33)  describes  the  convection  of  which  is  driven  by  a  source  term 
proportional  to  <f>i-  Eq.(34)  describes  the  reaction  of  i/q  on  $ j,  which  occurs 
via  a  modified  Poisson’s  equation.  In  the  limit  of  a  constant  density,  \ /q  is 
simply  proportional  to  the  linearized  charge  density,  and  Eqs. (33.34)  have 
a  clear  physical  interpretation,  the  transport  of  charge  and  its  subsequent 
modification  of  the  potential,  which  is  somewhat  obscured  by  the  presence 
of  the  density-gradient  dependent  terms. 

We  solve  Eqs. (33, 34)  for  the  unstable  modes  with  the  largest  growth- 
rates  by  integrating  numerically  from  random-noise  initial  conditions;  we 
then  measure  the  growth-rate  of  the  mode  which  emerges  out  of  the  noise 
The  numerical  method  is  straightforward;  Eq.(33)  is  solved  by  a  predictor- 
corrector  scheme,  with  4> i  obtained  from  Eq,(34)  through  a  tridiagonal  ma¬ 


trix  solution.  Throughout,  we  impose  the  conducting-wall  boundary  con¬ 
dition,  d>\{x  =  0)  =  0,  and  a  far-field  condition  <b\{x  =  L)  —  0,  where 
L  a,  d.  The  overall  method  is,  of  course,  limited  to  determining  the 
eigenvalues  of  the  unstable,  fastest-growing  modes;  but  this  is  all  we  re¬ 
quire  for  comparison  with  the  simulation  results. 

In  Fig. (10)  we  show  the  results  of  our  numerical  integration  of  Eqs.  (33,34 
with  the  parameters  of  Eqs. (26, 27).  We  have  plotted  the  the  real  and  imag¬ 
inary  parts  of  the  frequency  u  =  u>r  +  i 7  as  functions  of  kya.  From  Eq.(32) 
one  can  establish  the  general  functional  dependence  of  the  eigenfrequencies: 

u 1  =  —Sl(kya,kyd,cr)  (35) 

a 

where  Q  is  obtained  from  a  normalized  form  of  Eq.(32),  which  depends  only 
on  the  dimensionless  parameters  kya ,  kyd  and  a.  Furthermore,  the  roots 
occur  in  pairs  such  that  Q(  —  ky)  —  —fl‘(ky).  Thus,  for  a  given  choice  of 
d/a  and  <r,  one  needs  to  compute  numerically  the  dependence  of  Q  on  only 
one  free  parameter,  |&y|a,  and  for  only  positive  ky,  a  simplification  which 


Figure  10:  Growth-rates  and  frequencies  for  the  Kelvin-Helmholtz  instabil¬ 
ity,  with  the  equilibrium  profiles  of  Eqs. (21-25);  d/a  =  2,  a  =  0.73. 

reduces  the  numerical  work.  Fig. (10)  is  established  with  the  parameters  of 
Eqs.(26,27),  with  d/a  =  2  and  a  —  0.73. 

The  curve  for  the  growth  rate  7 (ky)  is  shown  in  Fig. (10a).  Its  shape 
is  roughly  parabolic,  with  cutoffs  at  kya  =  0  and  kya  =  0.9,  and  predicts 
a  maximum  growth  rate  of  7  =  O.OQOVo/a,  which  occurs  at  kya  =  0.43. 
The  curve  for  the  real  part  of  the  frequency  is  plotted  in  Fig. (10b).  This 
figure  shows  that  the  dispersion  relation  is  approximated  by  the  linear  re¬ 
lation  u>f{(ky)  =  —kyV0.  This  result  is  not  surprising:  the  unstable  Kelvin- 
Helmholtz  modes  can  be  thought  of  as  resulting  from  the  interaction  of 
two  counter-streaming  surface  waves,  moving  symmetrically  with  respect 
to  the  drift  velocity  at  the  inflection  point  in  the  velocity  profile[ll].  I11 
the  frame  co-moving  with  the  inflection  point,  the  interaction  leads  to  an 
almost  purely-growing  mode.  Now,  with  the  profile  of  Eq.(24),  the  velocity 
of  the  inflexion  point  is  precisely  —Vo,  and  this  results,  in  the  lab  frame, 
in  the  linear  dispersion  relation  noted  above.  The  fastest-growing  mode  is 
indeed  exactly  purely-growing  in  the  co-moving  frame,  and  results  in  the 
real  frequency  u:r  —  — 0.43V’o/a. 

In  Fig. (11)  we  plot  the  eigenfunction  Q\{x)  corresponding  to  the  max¬ 
imum  growth  rate.  In  this  plot,  we  can  see  that  while  the  peak  of  the 
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Figure  11:  Eigenmode  structure  for  the  fastest-growing  Kelvin-Helmholtz 
mode  of  Eq  (32) 

eigenftinction  occurs  at  x  —  1  5a,  that  is,  right  at  the  edge  of  the  shear 
layer,  the  eigenfunction  also  embodies  a  long  exponential  tail  which  pene¬ 
trates  deep  into  the  plasma,  up  to  a  distance  x  ~  10a. 

In  Figs. (12),  we  have  reproduced  Figs. (10),  but  have  normalized  all 
variables  to  the  physical  parameters.  Note  that  the  fluid  results  do  not  in¬ 
corporate  any  finite-gyro-radius  effects,  despite  what  aught  be  inferred  from 
the  presence  of  the  kyp,  term.  In  Figs. (12),  we  have  also  plotted  the  data  for 
the  growth  rates  and  the  real  parts  of  the  frequencies  as  directly  obtained 
from  Run  1.  The  comparison  shows  that  the  analytic  predictions  for  the 
real  part  of  the  frequency  agree  reasonably  well  with  the  simulation  results. 
On  the  other  hand,  while  the  growth-rates  for  the  long-wavelength  modes 
(m  =  1,2)  are  in  good  agreement,  fluid  and  simulation  results  diverge  as 
we  go  to  the  shorter  wavelengths  (m  =  3, 4,  5, 6):  while  the  fluid  model 
predicts  in  this  case  a  cutoff  at  kyp ,  =  1.4,  the  simulation  results  exhibit 
a  lower  cutoff,  at  kypi  —  0.9,  with  a  smaller  maximum  growth  rate.  Let 
7 F(ky)  denotes  the  growth-rate  predicted  by  the  fluid  theory,  and  7 y(ky) 
the  numerical  result  of  Run  1.  We  have  found  the  approximate,  empirical 
relation: 


7 N(ky)  as  ')f(ky)  -  0.1  (k2yp2)  (36) 

It  is  reasonable  to  think  that  the  “correction”  term  in  Eq.(36)  accounts  for 
finite-l armor  radius  effects  not  present  in  the  fluid  theory,  and  we  can  qual¬ 
itatively  account  for  the  quadratic  dependence  of  this  term  by  the  following 
derivation.  We  first  naively  add  the  orbit- averaged  contribution: 


A  2  t— '2 

VFLR  =  tPi  ^ 
4 


Ex: 
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to  the  ion  velocity  in  the  continuity  equation,  Eq.(13).  This  yields,  in  the 
limit  where  the  Navier-Stokes  equation  is  valid,  the  following  equation: 

Jj-VV  +  ^[4>, VV]  -  iu,clP|2W  =  0  (38) 

In  other  words,  there  is  now  a  diffusive  term  in  the  vorticity  equation.  This 
term  leads  to  an  additional  damping  v  which  has  the  dependence: 

v  ~  -udkl  +  fyp2  (39) 

and  we  note  that  the  quadratic  dependence  of  this  term  is  at  least  compat¬ 
ible  with  that  of  the  correction  term  in  Eq.(36). 
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Figure  12:  Comparison  of  the  analytic  results  for  the  Kelvin- Helmholtz 
instability  with  the  simualtion  results  of  Run  1;  (a)  growth  rates;  (b)  fre¬ 
quencies. 

6  Steady-State  Structure  of  the  Sheath 

In  the  previous  sections  we  discussed  the  transient  behavior  of  the  sheath; 
we  shall  now  examine  its  steady-state  structure.  We  shall  first  present  the 
numerical  results  of  Run  1,  and  discuss  the  structure  of  the  fields  and  of  the 
particle  densities.  We  shall  then  make  a  connection  with  the  fluid  theory 
of  Section  5,  by  fitting  a  solution  of  the  fluid  equations  to  the  observed 
potential  structures.  The  aim  of  this  procedure  is  to  make  explicit  the 
fluid  nature  of  the  vortices.  We  shall  close  this  section  with  a  discussion  of 
the  power  spectra,  which,  in  particular,  provide  information  on  the  short- 
wavelength  turbulence  which  complements  the  long-wavelength  vortices. 

6.1  Simulation  Results 

As  we  saw  in  Figs. (4),  after  u )ctt  ~  500  the  sheath  settles  into  a  steady-state, 
in  which  a  single  large  vortex  uniformly  drifts  along  the  wall.  This  large 
vortex  is  occasionally  accompanied  by  smaller  “satellites”,  which  reside  in 
the  system  for  intervals  of  order  uclT  ~  100.  The  average  drift  velocitv 
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uoy  of  the  principal  vortex  can  be  computed  from  the  slopes  of  the  tracks 
in  Fig. (4)  and  was  found  to  be  essentially  constant.  At  u>cti  =  700,  it  was 
found  to  be  |v0y|  /vti  =  0.44  (or  n0y  —  —0.065  in  the  units  of  ES2). 

A  three-dimensioned  perspective  plot  of  the  electrostatic  potential  <f>(x,y) 
at  u>dt  =  700  is  shown  in  Fig. (13).  This  figure  clearly  shows  that  the  vortex 
produces  a  very  sizeable  perturbation  of  the  edge  potentials.  In  Figs. (14) 
we  display  more  detailed  information  on  the  potential  of  Fig.(13),  in  the 
form  of  a  contour  plot  and  cross-sectional  plots.  The  maximum  poten¬ 
tial  drop  in  the  vortex  center,  relative  to  the  y-averaged  potential  profile 
(Fig. (14c)),  is  found  to  be  e8(j>/Ti  =  —2.4,  while  the  total  drop  of  the  y- 
averaged  potential,  from  the  wall  into  the  plasma,  is  eA 0/T,  =  —1.9.  The 
dimensions  of  the  separatrix  of  the  vortex  can  be  measured  from  Fig. (14a) 
and  are  found  to  be  ly  as  96  {ly/pi  ~  15),  and  lx  as  32  (lx/p,  %  5).  Thus,  if 
defined  by  its  separatrix,  the  vortex  has  the  shape  of  an  ellipse  elongated  in 
the  y-direction  (this  shape  is  distorted  by  the  choice  of  scales  in  Fig. (14a)). 
However,  the  funnel  of  the  vortex  flares  out  more  rapidly  in  the  y-direction, 
and  when  considered  say,  half-way  down  the  potential  well,  the  shape  of 
the  vortex  is  more  nearly  circular. 

In  Fig.(14d)  we  have  shown  the  y-averaged,  x-directed  electric  field 
Ex(x),  which  induces  the  drift  vy(x)  =  —Ex(x)/B0.  Fig.(14d)  shows  that 
the  y-averaged  velocity  profile  remains  strongly  sheared,  even  in  the  steady- 
state.  At  the  wall,  the  drift  velocity  is  ny(0)/e(t  =  1.25,  or  roughly  the 
ion  thermal  speed.  The  shear  layer  extends  to  about  three  ion  gyro-radii 
inward,  to  x  as  20. 

In  Fig. (15),  we  show  the  profiles  of  the  y-averaged  particle  densities. 
ne(x )  and  n,(x).  There  is  a  region  of  large  charge  separation,  in  which 
ne(x)  >  n,(x),  extending  to  x  as  20;  this  region  corresponds  to  the  electnc 
field  shown  in  Fig.(14d),  and  can  be  said  to  define  the  width  of  the  y- 
averaged  sheath.  Beyond  this  region,  for  20  <  x  <  32,  there  is  a  narrow 
quasineutral  region,  in  which  ne(x)  as  rc,(x). 

In  Figs. (16a, b),  we  have  displayed  the  “scatter  plots”  of  the  particle 
positions  at  u >cit  =  700.  A  striking  feature  of  these  scatter  plots  is  the 
existence  of  a  large  region  on  the  underside  of  the  vortex,  which  is  al¬ 
most  evacuated  of  particles.  We  shall  discuss  this  feature  in  Section  7. 
in  connection  with  our  discussion  of  particle  transport.  Cross-sections  of 
the  local  (unaveraged)  electron  and  ion  densities  at  u;c, t  =  700  are  shown 


Figure  14:  Details  of  the  potential  structure  of  Fig. (13),  u :c,t  =  700,  Run 
1;  (a)  contour  plot  of  <j)(x,y)\  (b)  y  cross-section  of  Fig. (a)  along  x  =  15, 
passing  through  the  vortex  center;  (c)  x  cross-section  of  Fig. (a)  along  y  = 
173,  passing  through  the  vortex  center;  (d)  profile  of  the  y-averaged  electric 
field  Ex(x). 


Figure  15:  Profiles  of  the  y-averaged  electron  and  ion  densities  at  u >ctt  = 
700,  Run  1. 


in  Figs. ( 16c, d),  onto  which  we  have  superposed  the  potential  profile  of 
Fig. (14b).  Let  us  note  the  presence  of  a  large  electron  density  at  the  core 
of  the  vortex  (Fig.(16c)).  As  can  be  seen,  there  is  no  corresponding  peak 
in  the  ion  density  (Fig.(16d)).  A  large,  net  negative  charge  at  the  center 
is  of  course  necessary  to  support  the  large  potential  drop  in  the  vortex: 
Figs. ( 16c, d)  provide  the  additional  information  that  this  charge  is  provided 
almost  entirely  by  trapped  electrons. 

In  Figs. (17)  we  look  at  the  structure  of  the  vortex  at  uclt  =  700, 
in  a  reference  frame  c^-moving  with  the  vortex,  at  a  (negative)  velocity 
vQy  =  —0.065.  The  potential  as  seen  in  this  reference  frame  is  found  by 
adding  to  the  potential  in  the  lab  frame  a  linear  potential  <t>c(x)  —  —  VyoBox, 
corresponding  to  the  electric  field  induced  by  the  uniform  motion  across 
the  background  magnetic  field.  Fig. (17c)  shows  that  Ex  has  an  approxi¬ 
mately  Unear  dependence  along  the  j-section  of  the  vortex,  with  maxima 
Ex/B0  «  ±1.5v*,.  The  field  component  Ey,  when  measured  along  the  y- 
section  of  the  vortex,  is  perhaps  better  approximated  by  a  cosine,  with 
maxima  Ey/Bo  zz  0.6vt,.  These  values  of  the  fields  give  us  an  estimate  of 
the  circulation  time  or  “bounce  time”  of  a  particle  around  the  potential 
weU  of  the  vortex.  Estimating  the  circumference  of  the  vortex  as  being 


Figure  16:  Distribution  of  particles  at  =  700,  Run  1;  (a)  electron 
scatter  plot;  (b)  ion  scatter  plot;  (c)  y  cross-section  of  the  electron  density 
ne{x,y )  along  x  =  15,  passing  through  the  vortex  center,  with  e<f>(x,y)/T, 
superimposed  (dashed  line);  (d)  identical  to  (c),  but  for  the  ion  density 
n,(x,y). 
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2  x  (/„/ 2)  «  15 pi  (we  are  estimating  the  circumference  half-way  down  the 
potential  well),  and  taking  the  mean  drift  velocity  to  be  v ^  ~  0.6etI,  we 
find  that  particles  will  have  a  bounce  period  rg  given  by  uiCiTg  %  25.  This 
indicates  that  the  bounce  frequency  ugx  of  the  electrons  is  of  order: 

^si  =  0.25  ujci  (40) 

where  the  estimate  is  for  electrons  roughly  half-way  up  the  potential  well 
of  the  vortex. 

6.2  Dependence  of  the  Sheath  Thickness  on  System 
Size 

In  Run  1,  the  vortex  grew  so  as  to  occupy  the  entire  width  in  x  of  the 
simulation  region,  and  the  nonneutral  sheath  also  extended  across  the  en¬ 
tire  region  (Fig.(15)).  Thus,  the  question  arises  whether  the  sheath  has 
a  “natural”,  self-limiting  width  in  x,  or  whether  it  can  grow  indefinitely, 
provided  the  system  is  large  enough  to  accomodate  it.  To  answer  this  ques¬ 
tion,  we  ran  a  simulation  in  which  Lx  was  2.5  times  larger  than  in  Run  1 
In  Figs. (18)  we  show  the  results  of  this  simulation,  otherwise  identical  to 
Rim  1,  but  with  Lx/ p,  =  13,  Ty/p,  =  20,  and  an  injection  rate  s  =  2.51. 
Figs. (18)  clearly  show  that  in  the  larger  system,  the  sheath  remains  local¬ 
ized  near  the  wall.  The  vortex  is  seen  to  extend  out  to  a  distance  lx  «  5 p,, 
beyond  which  the  fields  are  nonzero,  but  rapidly  decreasing.  In  Fig. (18c)  we 
have  shown  the  y- averaged  electric  field  Ex(x).  This  figure  shows  that  the 
shear  layer  is  more  localized  than  the  vortex  fields,  a  characteristic  which 
was  already  illustrated  in  Fig. (11),  where  we  saw  that  the  eigenfunctions  of 
the  Kelvin-Helmholtz  instability  extend  well  beyond  the  region  of  nonzero 
shear.  In  Fig.(18d)  we  show  the  profiles  of  the  steady-state  electron  and 
ion  densities.  This  plot  shows  that  in  agreement  with  Fig. (18c),  the  region 
of  charge  nonneutrality  is  confined  to  the  edge  layer,  0  <  x  <  5 p,.  Beyond 
this  edge  region  stretches  a  quasineutral  plasma,  which  has  a  density  profile 
that  is  roughly  parabolic. 
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Figure  17:  Cross-sections  of  the  fields  at  u>cif  =  700,  Run  1,  in  a  frame 
co-moving  with  the  vortex  at  a  velocity  vy  =  —0.44vu. 
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6.3  Localized  Solutions  for  the  Steady-State  Vortex 

The  characteristics  of  the  dominant  vortex  in  the  final  steady  state  were 
illustrated  in  Figs. (14).  We  Eire  now  faced  with  the  problem  of  providing  an 
analytic  solution  for  the  vortex  potential,  consistent  with  the  simulation  re¬ 
sults  and  with  the  cross-field  equations  (17,18).  Our  discussion  will  remain 
qualitative;  in  particular  we  shall  assume  that  in  our  situation  the  Navier- 
Stokes  equation  (21)  is  an  acceptable  approximation  to  the  full  cross-field 
equations  (17,18).  This  is  equivalent  to  assuming  that  the  effect  of  den¬ 
sity  variations  can  be  ignored,  by  taking  n  ss  constant.  While  there  are 
in  fact  strong  density  variations  in  the  sheath  (Fig.(16)),  we  are  adopting 
this  assumption  because  it  greatly  simplifies  our  quest  for  solutions  of  the 
cross-field  equations:  the  literature  on  the  Navier-Stokes  equation  is  large, 
and  many  analytical  solutions  are  available. 

The  solutions  of  the  Navier-Stokes  equation  which  are  available  in  the 
literature  are  essentially  of  two  types:  the  first  type  consists  of  local¬ 
ized,  soliton-like  solutions  (“modons”)[30];  the  second  of  periodic  solutions, 
formed  with  trains  of  nonlinear  waves.  We  shall  examine  these  in  turn. 

The  extent  of  the  vortex  of  Fig.(14a)  is  small  compared  the  system 
length  Ly ,  and  this  suggests  that  we  find  a  localized  solution  of  the  Navier- 
Stokes  equation  to  fit  to  the  observed  vortex  structure.  This  is  the  approach 
adopted  by  Horton  et  al.  [16],  who  derive  a  solution  for  the  asymptotic  state 
of  a  symmetric  shear  layer,  using  a  matching  technique  first  used  by  Sagdeev 
et  al.[25].  However,  we  have  not  had  much  success  with  this  approach:  this 
is  because  the  resulting  analytic  solutions  simply  do  not  reproduce  the 
large-amplitude  vortex  seen  in  our  simulations.  Below,  we  shall  first  review 
the  method  of  Sagdeev  et  al.  [25] .  We  shall  then  argue  that  a  more  general 
class  of  solutions  than  the  ones  used  by  Horton  et  al.  must  be  considered, 
if  we  are  to  approximate  the  steady-state  vortex  of  the  sheath.  We  shall 
finally  describe  a  simpler,  periodic  solution  (Stuart’s  solution),  which  at 
least  provides  a  rough  fit  to  the  observed  vortex  parameters,  leaving  to 
future  research  a  more  detailed  exploration  of  analytic  solutions  for  the 
asymptotic  states. 

We  now  review  the  general  procedure  for  obtaining  localized  solutions 
of  the  Navier-Stokes  equation[3l].  The  first  step  is  to  assume  a  structure 
stationary  in  a  reference  frame  moving  with  constant  velocity  vy  =  u.  To 
make  the  contact  with  Section  6.1,  we  write  u  =  — 1’0,  where  — 1’0  is  the 
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drift  velocity  of  the  vortex  as  seen  in  the  simulations.  Then,  defining  the 
transformed  potential  ip  through: 


<f>(x,y,t)  =  -v0x  +  ip(x,y  +  v0t)  (41) 

we  find  that  Eq.(21)  reduces  to  the  stationary  form: 

b,V20]  =  0  (42) 

where  the  Poisson  brackets  denote  the  operator  [/,  g]  =  dxfdyg  —  dyfdzg 
[27j.  To  solve  Eq.(42),  one  can  then  select  a  class  of  solutions  for  which: 

VV  =  /(«/-)  (43) 

where  f{ip)  is  a  given  function  of  U>.  If  anything,  Eq.(43)  suffers  from 
too  much  latitude  in  the  possible  choices  of  the  function  f(ip).  which  is 
essentially  arbitrary  in  the  form  given  above.  To  reduce  the  undeterminacy 
of  Eq.(43),  one  can  proceed  by  imposing  additional  constraints.  The  first 
simplification  is  to  assume  that  f(ip)  is  linear,  in  the  form: 

f{ip)  =  d  +  cip  (44) 

where  d  and  c  are  constants.  This  has  the  advantage  of  leading  to  an 
equation  which  is  analytically  solvable.  The  second  constraint  is  to  assume 
(43)  valid  only  in  a  restricted  region  of  the  plane,  say  in  a  circle  of  radius 
r  =  a,  and  require  that  the  solutions  of  (43)  match  to  a  prescribed  lami¬ 
nar  flow  outside  the  region.  Because  of  these  constraints,  many,  if  not  all 
of  the  arbitrary  constants  in  the  solution  of  the  interior  region  are  then 
determined. 

We  have  applied  the  method  outlined  in  the  preceeding  paragraph  in  two 
ways.  In  the  first,  we  attempted  to  model  the  edge  region  by  a  monopolar 
vortex,  symmetric  about  the  inflection  point  of  the  shear  layer.  In  the 
second,  we  attempted  to  fit  a  dipolar  vortex[32]  to  the  edge  region,  with 
the  wall  as  symmetry  plane  between  the  two  vortices  in  the  dipole,  and 
with  one  of  the  vortices  a  virtual  image  of  the  other  one,  resident  in  the 
plasma.  In  each  case,  the  interior  region  of  the  vortex  was  modelled  by 
Eq.(44)  with  d  =  0  and  c  =  —  k2,  resulting  in  the  Helmholtz  equation: 


V20  +  k2ip  =  0 


(45) 


The  matching  condition  at  the  boundary  between  the  vortex  interior  and 
the  unperturbed  exterior  flow  then  results  in  a  “modon  dispersion  relation" 
for  the  internal  wavevector  k. 

These  approaches  unfortunately  failed,  in  that  the  solutions  predicted 
vortices  with  small  amplitudes,  e\8<j)\  /T,  <C  1,  in  contradiction  with  the 
observed  e  |5<>|  /T,  ~  2.  The  weakness  of  the  vortices  obtained  in  the  ana¬ 
lytic  solutions  was  linked  to  our  choice  of  the  smallest  root  of  the  dispersion 
relation  for  k,  this  choice  leading  to  a  flow  which  closely  approximates  the 
external,  unperturbed  flow.  By  choosing  the  next,  larger  root,  one  does  ob¬ 
tain  a  large-amplitude  vortex,  with  e  |£d>|  /T,  ~  0(1).  However,  this  choice 
also  leads  to  strongly  reversed  flows[l6],  which  do  not  ressemble  the  circular 
flows  seen  in  the  simulation  vortices. 

We  suspect  that  the  failure  of  the  method  is  tied  to  the  choice  of  a 
linear  function  for  because  such  a  choice  does  not  produce  solutions 

which  are  sufficiently  “self- binding”.  As  an  alternative,  one  might  consider 
an  exponential  dependence  in  the  function  f(tp),  with  the  ansatz: 

f(xp)  ~  aexp(60)  (4G) 

For  values  of  ip  of  order  unity,  the  exponential  can  provide  the  necessary 
curvature  in  the  solution  to  ip  so  that  within  a  localized  region,  il>  can  match 
to  some  external  flow. 

We  have  not  attempted  to  solve  Eq.(43)  with  the  exponential  depen¬ 
dence  /( te)  ~  aexp(btp)  for  a  strictly  localized  solution.  Rather,  we  have 
fallen  back  onto  a  periodic  solution  existing  in  the  literature,  which  does 
satisfy  Poisson’s  equation  with  this  form  of  This  is  Stuart's  solu¬ 

tion  .  which  can  be  written: 

ib  =  —  log  [cosh(£.r)  +  .4cos(£j/)]  (47) 

K 

In  Eq.(47),  the  wavenumber  k  is  adjustable  and  determines  the  periodicity 
of  the  solution,  with  wavelength  A  =  2w/k.  The  parameter  A  is  also  ad¬ 
justable,  with  A  <  1,  and  determines  the  maximum  depth  of  the  potential 
well.  One  can  verify  that  the  function  ip  of  Eq.(47)  satisfies: 


V2ip  =  kv0  exp(  —  2k\b/vo)  (48) 

This  corresponds  to  a  charge  density  with  a  strong  dependence  on  the 
potential: 


p(rp)  =  —kv0  exp(—2kip  /  v0)  (49) 

Note  that  if  we  normalize  x  and  ip  in  Eq.(48)  according  to  £  =  kx  and 
£  =  kip/ Vo,  then  Eq.(48)  reduces  to  the  simple,  general  form: 

V|C  =  exp(-2C)  (50) 

While  Stuart’s  solution  is  an  exact  stationary  solution  of  the  Navier- 
Stokes  equation,  its  choice  is  a  rough  physical  guess  insofar  as  it  is  not 
the  exact  consequence  of  the  time-dependent  problem.  Keeping  this  spirit 
of  approximation  in  mind,  we  shall  now  determine  the  free  parameters  in 
Stuart’s  solution  by  matching  to  the  observed  simulation  results. 

6.4  Periodic  Solution  for  the  Steady-State  Vortex 

Using  Eq.(41),  we  transform  the  potential  ip  of  Eq.(47)  back  to  the  potential 
(p  in  the  stationary  laboratory  frame: 

<P{x,y)  =  - VqX  +  j^\og(cosh(k(x  -  b ))  +  Acos(k(y  -  v0t)))  (51) 

In  this  equation,  we  have  also  chosen  a  new  line  of  symmetry  in  x ,  by 
shifting  the  vortex  centers  to  x  =  b.  Note  that  Eq.(51)  predicts  vy(x )  — ♦  0 
for  k(x  —  b)  1,  and  very  roughly  vy{x  =  0)  —  2v0  at  the  wall.  By 

inspection  of  Fig. (14),  we  find  that  the  vortex  center  lies  at  b  =  16  (6/p,  = 
2.5).  For  the  wavelength  k,  we  have  k  =  2n/Ly  =  0.0245  (kp,  =  0.155). 
The  vortex  velocity  v0  was  already  noted  in  Section  6.1,  where  it  was  fomid 
to  be  volvn  =  0.44. 

To  determine  the  remaining  adjustable  parameter  .4,  we  consider  the 
maximum  depth  of  the  potential  well  along  the  line  of  symmetry,  x  =  b.  In 
physical  units,  this  is: 


or,  solving  Eq.(52): 


With  the  values  derived  above  and  e  \8<t>\  /Tt  v  2.4,  we  find  A  %  0.57. 
To  summarize,  we  have  found  that  Eq.(51)  gives  a  rough  estimate  for  the 
steady-state  vortex  with  the  parameters: 


vo/vtl  =  0.44  kpi  =  0.155  (54) 

6/p,  =  2.5  A  =  0.57  (55) 

Let  us  examine  some  properties  of  Stuart’s  solution  with  these  param¬ 
eters.  With  k  |z  —  6|  1,  Eq.(51)  predicts  the  potential: 

4>(x,y)  =  -v0x  +  y  log(cosh(fc(x  -  6)))  (56) 

and  hence  the  flow: 


vy(x)  =  — <f>(x,y )  =  t’o  (tanh(fc(x  -  6)) 


1) 


(57) 


The  analytic  form  of  this  flow  is  identical  to  that  which  we  assumed  for  the 
initial  conditions  of  the  linear  stability  analysis  of  Section  5.2: 


«oy(x)  =  V0(tanh(x/a  -  1)  -  1)  (58) 

provided  we  take  Vo  =  u0.  However,  the  final  value  of  k  observed  in  the 
simulation  imposes  a  much  broader  scale  length  than  the  initial  gradient 
length  a,  as  we  have  1/k  =  6.5p,,  which  is  much  larger  than  a  =  0.4p;. 
From  Eq.(51)  we  can  obtain  the  total  width  in  x  of  the  separatrix.  This  is 
given  by: 

h  =  ^-cosh-1(l  4-  2,4)  (59) 

We  find  that  h  ss  8p,,  which  is  somewhat  larger  than  the  observed  vortex 
diameter  lz  5p, .  If  we  find  this  discrepancy  disturbing,  we  should  perhaps 
instead  fink  the  vortex  size  not  to  the  extent  of  its  separatrix,  but  to  the 
scale-length  of  its  (/-averaged  flow,  which  is  somewhat  sharper.  Averaging 
over  the  y  coordinate  at  any  fixed  x,  we  find  that: 


vy(x)  =  -v0  -  Ko 


sinh(fc(x  —  6)) 


(cosh2(fc(x  —  6))  —  A2) 


1/2 


(60) 
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With  A  =  0.57,  we  find  that  vy(x)  achieves  3/4  of  its  maximum  value  at  a 
distance  /  =  1.3p,  from  the  plane  of  symmetry.  We  shall  now  decree  that 
this  is  the  new  scale-length  for  the  shear  layer,  rnd  that  the  line  of  symmetry 
has  shifted  in  such  a  way  as  to  accomodate  it,  by  defining  b  =  l  ss  1.3p,. 
For  the  parameters  of  Simulation  1,  this  yields  b  ss  8.1,  which  is  close  to 
the  observed  value  of  the  coordinate  of  the  inflection  point. 

The  value  of  A  which  we  found  above,  A  =  0.57,  is  comparable  to 
that  found  when  one  fits  the  model  of  Eq.(51)  to  the  simulation  results  of 
Corcos[20,21,22],  In  a  simulation  of  the  Kelvin-Helmholtz  instability,  with 
a  tanh  profile  for  the  shear  layer  similar  to  the  one  above,  Corcos  initially 
excited  the  mode  with  ka  =  0.45.  The  saturated  vortex  which  developped 
from  the  initial  perturbation  was  seen  to  have  a  separatrix  width  satisfying 
kh  =  10.0,  which  from  Eq.(59)  corresponds  to  ,4  =  0.45 

Finally,  to  measure  the  vorticity  extracted  by  the  instability  from  the 
shear  layer,  we  consider  the  circulation  T  =  J  uodxdy.  The  ratio  of  the  total 
circulation  in  the  model  vortex  of  Eq.(51)  to  the  initial  circulation  T0  in  a 
wavelength  2it/k,  is  given  by: 

=  —  tan-1  Va  =  0.82  (61) 

T0  n 

so  that  the  simulation  results,  and  the  analytic  model  for  the  vortices, 
are  compatible  with  a  situation  in  which  most  of  the  vorticity  (and  hence 
most  of  the  net  charge  density  in  the  shear  layer)  resides  within  the  vortex 
separatrix. 

6.5  Power  Spectra 

We  shall  now  discuss  the  frequency  and  wavenumber  dependence  of  the 
power  spectra.  As  in  Section  4.3  we  consider  the  transforms  of  the  potential 
<p(x,y,t)  at  a  fixed  x,  located  in  the  midplane  of  the  simulation  region. 
x  =  Lx/ 2,  but  we  now  transform  in  time  as  well  as  in  y.  More  precisely,  for 
each  wavenumber  ky,  the  power  spectrum  P{ky,  u>)  is  defined  as  the  Fourier 
transform  in  time  of  the  autocorrelation  function  of  the  mode: 


P(ky,  u)  =  [  dr  e,UJT  <  <j>(ky,t)d>m(ky,t  +  t)  > 


(62) 


In  Eq.(62)  the  average  <  . . .  >  is  obtained  by  performing  a  time-average  on 
the  simulation  data  from  the  time-interval  0  <  <  1075.  The  resulting 

autocorrelation  function  is  also  passed  through  a  lag  window [33],  resulting 
in  a  smoothed  spectrum. 

Let  us  first  consider  the  total  spectral  energy  for  each  Fourier  mode 
4>(ky,t).  This  is  simply: 

<  \<t>(ky,t)\2  >=  —  P(kv,w)  (63) 

In  Fig. (19)  we  plot  <  \4>(ky,  f)|2  >  and  compare  its  values  to  <  |ri(ky,f)|2  >th. 
the  distribution  of  energy  in  thermal  equilibrium.  The  theoretical  expres¬ 
sion  for  the  spectrum  in  thermal  equilibrium  is  given  by: 

<  \i(k„,t)\*  >“=  ~jh.  (>  -  (1  +  (W!r1,!)  (64) 

which  is  obtained  from  the  fluctuation-dissipation  theorem,  assuming  an 
infinite,  homogeneous  plasma[34].  While  Eq.(64)  does  not  account  for  all 
numerical  discreteness  effects,  nor  for  the  boundedness  of  the  system,  we 
believe  that  it  provides  a  qualitative  estimate  for  the  level  of  thermal  fluctu¬ 
ations.  The  first  notable  feature  of  Fig. (19)  is  that  the  spectrum  is  peaked 
at  low  wavenumbers,  and  is  strongly  cutoff  beyond  kyp,  =  1,  a  feature 
already  observed  in  Fig. (7),  where  we  showed  a  snapshot  of  modal  ampli¬ 
tudes.  The  second  distinctive  feature  of  Fig. (19)  is  that  the  values  of  the 
power  spectrum  Eire  quite  larger  than  those  of  the  thermal  spectrum  over  a 
large  range  of  wavenumbers,  up  to  kyp,  «  3.  This  latter  feature  indicates 
that  despite  their  smaller  amplitudes  (as  compared  to  the  ed>/T,  ~  1  vor¬ 
tices),  the  short  wavelength-modes  are  also  collective  plasma  oscillations, 
and  cannot  be  assimilated  with  the  thermal  fluctuations. 

We  shall  now  consider  the  power  spectra  in  more  detail,  as  a  function  of 
frequency  u>.  In  Figs. (20)  we  have  plotted  the  spectra  P(ky,u>)  for  the  first 
7  wavenumbers  (omitting  ky  =  0),  covering  the  range  0.155  <  kypi  <  1.085. 
A  first  feature  of  these  figures  is  that  in  each  graph,  the  spectrum  peaks  at 
a  different  value  of  ui,  w  =  u ip,  suggesting  a  dispersion  relation  of  the  form 
uiP  =  uj(ky).  Now,  if  the  modes  are  simply  static  structures,  carried  along 
at  the  vortex  speed  vy  —  —  Vo,  then  we  should  have  the  dispersion  relation: 


Figure  19:  Power  spectrum  for  Rim  1,  time-averaged  over  0  <  uc,t  <  1500 
(full  line).  We  have  also  plotted  the  estimated  thermal  spectrum  (dashed 
hne).  All  data  is  taken  along  the  y  cross-section  at  x  =  Lxj 2  =  16. 

In  Fig. (21),  we  compare  this  dispersion  relation  to  the  values  of  u Jp{ky):  we 
can  see  that  the  agreement  is  reasonably  good.  Let  us  note  a  second,  quali¬ 
tative  feature  of  Fig. (21),  which  we  believe  reveals  a  fundamental  aspect  of 
the  edge  turbulence:  this  is  the  progressive  broadening  of  the  spectra  as  we 
move  to  larger  kypt.  While  the  modes  with  kyp,  <C  1  have  a  spectrum  which 
is  both  narrow  and  confined  to  |u>|  uci,  the  modes  with  kyp,  ~  1  have  a 
broad  spectrum,  with  peak  near  u>cl  and  with  width  |Au>|  ~  u>c,.  There  is 
also  more  structure  in  the  short-wavelength  modes,  with  the  amplitude  of 
many  sidebands  comparable  to  that  of  the  main  peak. 

These  qualitative  observations  on  the  structure  of  the  spectra  suggest 
that  the  edge  fields  can  be  partitioned  into  two  physically  distinct  compo¬ 
nents,  which  can  be  described  as  follows: 

1.  There  is  a  “coherent”  component  of  the  fields,  with  \ky\p,  <C  1  and 
|w |  <  u ;CI,  corresponding  to  the  large  and  stable  structures  seen  in 
the  plasma  -  the  vortices.  This  component  arises  out  of  the  inverse 
cascade  of  wavenumbers  which  occurs  during  the  transient  build-up 
of  the  svstem. 
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Figure  20:  Power  spectra  in  Run  1,  time-averaged  over  0  <  u :c,t  <  1500. 
for  Fourier  modes  m  =  1  to  4;  all  data  is  taken  along  the  y  cross-section  at 
x  =  Lx/2  =  16  ( continued  on  the  next  page). 
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Figure  20:  ( continued  from  the  previous  page )  Power  spectra  in  Run  1. 
time-averaged  over  0  <  u ;cit  <  1500,  for  Fourier  modes  m  =  5  to  7:  all  data 
is  taken  along  the  y  cross-section  at  x  =  Lx/ 2  =  16. 
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Figure  21:  Dispersion  relation  obtained  from  Figs. (20);  the  dashed  line 
denotes  the  Unear  dispersion  relation  |u>|  =  ky  |t?0|  where  v0  =  —0.44vt,  is 
the  vortex  velocity. 

2.  There  is  also  a  “turbulent”  component  of  the  fields,  with  \ky\ p,  ~  1 
and  |u>|  ~  wCj,  corresponding  to  the  small,  short-Uved  fluctuations. 
This  component  arises  out  of  a  forward  cascade  of  wavenumbers, 
which  we  shall  qualitatively  discuss  in  the  next  section. 

The  partition  scheme  outhned  above  is  based  on  qualitative  considera¬ 
tions  which  we  have  adopted  in  the  absence,  on  our  part,  of  a  more  precise 
theoretical  understanding.  In  Section  7.5,  we  shall  nonetheless  use  this 
scheme  to  estimate  diffusion  coefficients,  thereby  obtaining  at  least  semi- 
quantitative  results.  We  leave  to  future  research  the  development  of  a  more 
quantitative  understanding  of  the  mechanisms  underlying  the  edge  turbu¬ 
lence. 

6.6  Source  of  the  Short-Wavelength  Fluctuations 

What  are  the  mechanisms  underlying  the  partition  of  the  edge  turbulence 
into  a  spectrum  with  a  mixture  of  coherent  and  turbulent  components?  The 
answer  is  related  to  that  of  general  problems  in  two-dimensional  turbulence. 


in  particular  to  the  organization  of  forward  and  inverse  cascades  of  energy 
in  wave-number  space  Thus,  we  fear  that  a  quantitative  answer  to  the 
question  probably  requires  a  range  of  investigation  largely  outside  the  scope 
of  the  present  paper.  However,  in  keeping  with  the  spirit  of  the  analysis  of 
the  previous  sections,  we  shall  at  least  suggest  some  mechanisms  by  which 
the  short  wavelength,  “turbulent”  components  can  arise  in  the  system. 

Our  point  of  view  is  that  the  quasi-stable  vortex  state  observed  in  the 
simulations  can  be  regarded  in  a  first  approximation  as  a  completely  stable 
steady-state  of  the  form  given  by  Eq.(51).  The  problem  is  then  to  predict 
how  short-wavelength  ( kypj  ~  1),  high-frequency  (u >  ~  u>ct )  fluctuations  can 
arise  from  such  a  nonlinear  equilibrium  state.  The  analysis  is  complicated 
by  at  least  two  general  factors.  The  first  is  that  the  equilibrium  state  is 
itself  a  nontrivial,  nonlinear  solution  to  the  cross  field  equations,  and,  for 
instance,  linearizing  about  this  state  is  a  procedure  greatly  complicated 
by  the  specific  vortex  geometry.  The  second  factor  of  complication  arises 
because  we  are  considering  a  regime  for  which  the  fluid  equations  are  no 
longer  valid,  and  thus  the  full  kinetic  description  is  presumably  necessary. 
At  the  very  least,  an  equation  of  the  form  of  Eq.(38)  (but  with  a  rigorous 
derivation  of  FLR  effects  to  order  ( kp ,)2)  will  have  to  be  considered. 

Within  the  general  context  just  outlined  above,  we  shall  be  content  to 
suggest  the  following  mechanisms  for  the  excitation  of  the  short-wavelength 
fluctuations: 

1.  Instabilities  in  the  shear  layer  outside  of  the  vortex.  As  we  saw  in 
Figs. (4),  every  now  and  then  satellite  vortices"  are  generated  at 
some  distance  from  the  main  vortex  and  subsist  for  a  time  period 
ujCit  ~  100.  This  suggests  a  continuous  generation  of  small-wavelength 
modes  at  some  distance  from  the  main  vortex;  presumably  these 
modes  are  then  convected  up  and  down  the  shear  layer  relative  to 
the  main  vortex,  and  eventually  merge  with  it. 

2.  Instabilities  inside  the  main  vortex:  In  viscous  hydrodynamics,  large 
eddies  are  rarely  absolutely  stable,  and  are  seen  to  generate  smaller 
eddies  through  viscous  shear.  This  suggests  that  we  investigate  the 
stability  of  the  vortex  to  radial  eigenmodes  with  angular  dependences 
of  the  form  e,md .  As  we  shall  see  in  Section  7.1,  the  bounce  frequency 
of  electrons  at  the  bottom  of  the  potential  well  of  the  vortex  is  of  order 
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u>b  ~  0.5u)ct.  The  excitation  of  modes  at  harmonics  of  this  frequency 
could  account  for  the  presence  of  the  high  frequency  fluctuations. 

3.  Trapped-electron  instabilities:  The  presence  of  trapped  electrons  in 
the  vortices  suggests  the  possibility  of  sideband  instabilities  as  studied 
by  Kruer  and  Dawson[35]  in  connection  with  large-amplitude  plasma 
waves. 

4.  General  parametric  instabilities:  Rather  than  consider  the  stability 
of  a  linearization  around  the  vortex  equilibrium,  we  can  consider  the 
more  general  question  of  the  stability  of  arbitrary  three-wave  pro¬ 
cesses.  In  particular,  there  is  the  possibility  of  three-wave  coupling 
to  the  ky  =  0  mode,  a  process  analysed  in  connection  with  drift 


7  Transport 

An  essential  feature  of  the  dynamics  of  the  cross-field  sheath  is  the  exis¬ 
tence  of  an  outward,  ambipolar  transport  of  particles.  In  the  steady-state, 
this  transport  is  nearly  constant  in  time,  and  equal  to  the  net  influx  due 
to  pair  creation.  What  is  the  driving  mechanism  of  this  transport?  In  the 
previous  section,  we  saw  that  the  power  spectrum  could  be  roughly  par¬ 
titioned  into  two  qualitatively  different  contributions:  a  long-wavelength, 
“coherent”  component,  which  is  associated  with  the  long-lived  vortices;  and 
a  short-wavelength  “turbulent”  component,  associated  with  the  short-lived 
fluctuations  which  satisfy  the  condition  |fcy|pt-  ~  1.  To  this  spectrum,  we 
must  add  an  overall  thermal  background,  due  to  the  excitation  of  waves  by 
the  discrete  particles.  In  the  present  section,  we  shall  estimate  the  contri¬ 
bution  to  transport  from  each  of  these  components. 

The  “coherent”  component  of  the  fields  is  composed  of  the  large,  long- 
lived  vortices.  Its  representation  is  that  of  a  uniformly  drifting,  constant- 
amplitude  vortex  of  the  form: 

<f>co/*(x,y,f)  =  $(x,y  -  ut)  (66) 

where  u  =  —  v0  is  the  vortex  drift-velocity.  Let  us  consider  the  effect  of 
(66)  on  transport.  First,  because  of  the  functional  dependence  of  <I\  we 
can  transform  to  a  reference  frame  co-moving  with  the  vortex,  in  which  the 
potential  becomes  strictly  static.  Thus,  we  need  only  consider  the  effect  of 
the  time-independent  potential  $(x,y). 

Let  us  first  consider  the  effect  of  $>(x,y)  on  the  ion  motion.  On  account 
of  their  large  gyro-radii,  the  ions  would  appear  to  be  a  prion  the  least 
well-confined  species.  Furthermore,  a  potential  such  as  (66)  can  in  princi¬ 
ple  create  outward  transport,  by  moving  the  ions  in  intrinsically  stochastic 
orbits.  This  judgment  is  based  on  the  observation  that  the  ion  gyroradius 
is  an  appreciable  fraction  of  the  vortex  size  (  p,//x  ~  1/5  ).  suggesting 
that  the  adiabaticity  of  the  ion  E  x  B  drifts  is  destroyed  in  the  presence 
of  the  vortices,  with  stochastic  motion  as  a  result.  Of  course,  a  quantita¬ 
tive  analysis  of  such  a  transport  mechanism  requires  rigorous  Hamiltonian 
methods,  as  for  instance  in  [37].  However,  we  would  like  to  argue  that 
such  an  analysis  of  ion  transport  is  not  immediately  necessary,  because, 
without  a  ronoommitant  outward  electron  flow,  whatever  outward  ion  flow 


might  occur  through  intrinsic  stochasticity  is  strictly  self- limiting.  This  is 
because,  as  the  ions  are  driven  to  the  wall,  they  build-up  a  positive  charge 
on  its  surface,  which  eventually  forms  an  impassable  potential  barrier.  It 
is  not  hard  to  estimate  that  the  required  electric  field  for  stopping  almost 
all  ions  is  of  order  Ex/ Bo  ~  vti  (see  Eq.(  104)). 

Now,  for  the  electrons,  there  is  almost  no  outward  transport  in  the 
presence  of  the  potential  $(x,y).  If  $(x,y)  is  rigorously  time-independent, 
then  the  electron  orbits  will  remain  regular  (because  pe  <C  hortcx)-  Through 
their  Ex  B  motion  the  electrons  will  closely  follow  the  equipotentials  $(x,  y) 
at  all  times,  and  the  great  majority  of  electrons  will  remain  forever  confined, 
by  being  locked  onto  a  given  equipotential.  The  exception  will  be  for  those 
electrons  on  equipotentials  which  come  close  to  the  wall,  within  a  layer  of 
width  pe,  within  which  the  finite  electron  gyro-radius  enables  them  to  strike 
the  wall.  However,  this  class  of  equipotentials  defines  a  “scrape-off”  layer 
which  is  still  very  narrow  compared  to  the  system  size,  and  whose  existence 
does  not  explain  the  outward  transport  of  electrons  which  are  at  a  distance 
x  >  pe  from  the  wall. 

Thus,  we  must  conclude  that  by  itself,  the  “coherent”  component  of  the 
spectrum  is  unable  to  induce  transport.  Furthermore,  the  arguments  given 
above  suggest  that  this  situation  is  primarily  due  to  the  lack  of  electron  mo¬ 
bility.  While  it  appears  that  the  ions  might  have  an  intrinsic  diffusion  tied 
to  their  stochastic  orbits,  the  requirement  of  ambipolarity  automatically 
suppresses  their  ou‘  ard  transport,  whenever  the  electrons  themselves  re¬ 
mained  confined.  Thus,  the  electrons  form  a  “bottleneck"  in  transport.  In 
what  follows,  we  shall  therefore  concentrate  on  the  electron  diffusion  as  the 
dominant  transport  mechanism,  inducing  loss  of  particles  of  both  species: 
we  shall  assume  that  the  ions  easily  “follow  the  electrons",  and  that  their 
outward  flux  is  simply  equal  to  the  independently  determined  electron  flux. 

The  failure  of  the  coherent  part  of  the  spectrum  to  induce  transport 
suggests  that  the  turbulent  component  of  the  spectrum  plays  the  central 
role,  by  creating  an  effective  collisionality,  resulting  in  a  diffusion.  Because 
of  the  short- wavelength  nature  of  the  fluctuations,  we  suspect  that  this 
diffusion  coefficient  is  of  a  local  nature,  resulting  from  scattering  events 
over  a  length  scale  of  order  p,  or  smaller.  However,  we  should  note  that  in 
such  a  situation,  the  coherent  spectrum  might  still  have  a  sizeable  effect. 
This  consideration  arises  because  though  the  coherent  spectrum  does  not 
induce  transport  by  itself,  it  might  enhance  this  diffusion,  by  behaving  as 


a  convection  cell.  In  this  latter  scenario,  transport  is  enhanced  because  the 
vortex  can  ferry  particles  from  the  interior  to  regions  much  closer  to  the 
wall,  on  a  time-scale  shorter  than  the  global  diffusion  time  in  the  presence 
of  only  the  fluctuations. 

An  important  question  to  be  first  decided  is  whether  the  dominant 
source  of  transport  is  due  to  the  turbulent  collective  modes,  or  whether 
it  is  a  consequence  of  collisional  effects.  We  shall  consider  this  question  in 
Section  7.4,  where  we  show  that  the  scaling  of  transport  in  the  system  is  not 
collisional.  However,  we  shall  first  piesent  some  numerical  measurements 
of  transport  as  obtained  in  Run  1. 

7.1  Orbits  of  Test  Particles 

In  ordei  to  gauge  the  rate  of  transport  of  particles  across  the  system,  we 
set  up  a  numerical  experiment  in  which  we  followed  the  motion  of  a  set 
of  test  particles,  composed  of  five  electrons  and  five  ions  which  were  ini¬ 
tially  placed  at  the  bottom  of  the  vortex.  In  Fig. (22a).  we  show  the  initial 
particle  positions  at  u >cit  =  1075.  We  followed  the  motion  of  the  particles 
over  the  time  interval  1075  <  ~CIf  <  1575,  the  final  positions  being  shown 
in  Fig. (22b).  The  test  particles  moved  in  response  to  the  electric  fields  in 
the  simulation  region,  but  were  completely  passive,  in  that  they  did  not  in 
turn  affect  the  fields  through  Poisson’s  equation.  Thus,  their  presence  did 
not  modify  the  dynamics  of  the  system.  Furthermore,  in  order  to  follow 
as  closely  as  possible  the  guiding-center  motion  of  the  electrons,  we  intro¬ 
duced  an  unphysical  damping  in  the  mover  of  the  test  electrons[38].  which 
suppressed  gyro-motion,  and  retained  only  the  guiding-center  motion.  As 
we  shall  see,  a  consequence  of  this  modification  of  the  mover  is  that  the 
test  electrons  remain  in  the  system  much  longer  than  the  field  electrons, 
because  they  have  to  move  closer  to  the  wall  to  be  scraped-off.  A  rough 
indication  of  the  better  confinement  of  the  test  electrons  can  be  seen  in 
Fig. (22b),  which  shows  that  while  by  u jc,t  =  1575  most  of  the  ions  have 
impacted  into  the  wall,  most  of  the  electrons  are  still  confined.  Without 
the  damping  of  their  gyro-motions,  the  test  electrons  are  lost  at  the  same 
rate  as  the  test  ions. 

In  Figs. (23),  we  show  the  time  histories  of  the  x  coordinates  of  the 
five  test  electrons,  with  a  blow-up  of  the  motion  of  Electron  1  shown  in 
Fig. (24).  The  plots  show  that  all  electrons  are  initially  trapped  at  the 
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a)OJc(t-1075 


b)CUcjt-1575 


o  :  Electrons 
a  :  Ions 


Figure  22:  Positions  at  two  different  times  of  a  set  of  test  particles,  five 
electrons  and  five  ions,  in  Run  1.  These  are  superposed  on  the  correspond¬ 
ing  equipotential  contours,  (a)  u :cit  =  1075.  initial  positions  of  the  particles, 
with  electrons  and  ions  superposed;  (b)  u,'c,t  =  1575:  final  positions,  includ¬ 
ing  the  particles  which  have  impacted  into  the  wall. 
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Figure  24:  Enlargement  of  the  time-history  of  xe(t)  for  electron  1  of 
Fig. (23);  1075  <  uclt  <  1325,  Run  1. 

bottom  of  the  vortex,  and  remain  trapped  for  a  considerable  number  of 
bounce  oscillations.  From  Fig. (24),  we  can  estimate  the  bounce  frequency 
of  the  electrons  near  the  bottom  of  the  vortex;  it  is  found  to  be: 

uB2  =  0.5^  (67) 

This  estimate  is  to  be  compared  with  that  of  Eq.(40),  u>bi  =  0.25u>CI,  in 
which  we  evaluated  u >b  from  the  measured  values  of  the  E  x  B  drift  inside 
the  vortex.  As  ljbi  was  evaluated  for  a  particle  roughly  half-way  up  the 
potential  well,  it  is  not  surprising  that  u>B2  >  ^bi,  indicating  faster  bounce 
oscillations  at  the  very  bottom  of  the  vortex. 

Figs. (23)  show  that  over  many  bounce  oscillations,  the  electrons  be¬ 
come  progressively  untrapped,  by  migrating  toward  the  sepaxatrix  of  the 
vortex.  Thus,  xe(t)  exhibits  larger  vertical  excursions  as  time  progresses, 
and  the  bounce  oscillations  slow  down  towards  the  end  of  the  duration  of 
the  trapped  orbits.  After  reaching  the  separatrix,  the  electrons  become  un¬ 
trapped,  and  then  undergo  several  kinds  of  motion:  for  instance,  Electron 
1  begins  to  execute  large  excursions  across  the  system,  over  a  time-scale 
ucit  ~  100;  Electron  2  is  almost  immediately  driven  to  the  wall;  Electron  5 


Figure  25:  Parametric  plots  of  (xe(t),  ye(t))  for  electrons  p  =  1,2  of 
Figs. (23);  1075  <  uclt  <  1325,  Run  1. 


first  behaves  as  Electron  1  and  displays  large  excursions  across  the  system, 
but  is  eventually  retrapped  by  the  vortex  at  u>c,t  as  1450.  Note  that  in 
almost  all  cases,  the  electrons  are  driven  within  a  distance  pe  of  the  wall, 
and  without  the  damping  of  their  gyro-motion  would  have  been  scraped-off 
from  the  plasma. 

In  Figs. (25)  we  show  parametric  plots  of  the  electron  position  (xe(f),  ye(t)). 
so  as  to  give  a  general  feeling  for  the  electron  motion.  In  these  figures,  one 
can  distinguish  the  trapped,  oscillating  orbits,  from  the  relatively  straight, 
untrapped  trajectories. 

In  Figs. (26),  we  show  the  motion  of  test  Ion  1,  by  plotting  the  time 
histories  of  x,(t),  y,{t)  and  (*,-(<),  y<(f)).  Though  the  ion  is  initially  trapped 
in  the  vortex,  as  can  be  seen  in  Fig.(26a),  its  motion  when  trapped  is  less 
distinctive  than  for  the  electrons.  Thus,  in  Fig.(  2(,c)  we  do  not  see  the  great 
loops  that  were  seen  in  Figs. (24).  This  is  because  the  large  gyro-orbit  of  the 
ion  blurs  the  motion  of  the  guiding-center,  and  also,  presumably  destroys 
its  exact  E  x  B  drift.  A  distinctive  feature  of  Fig. (26c)  is  the  repeated 
repulsion  of  the  ion  by  the  strong  electric  field  in  the  sheath;  this  is  seen  as 
a  set  of  almost  point-like  reflections  of  the  ion  orbit. 

While  Figs. (23)  will  provide  us  in  Section  7.3  with  an  estimate  of  the 

70 


.-v- 


average  global  diffusion  time,  they  do  not  indicate  in  which  specific  regions 
transport  is  greatest.  We  can,  however,  obtain  an  indication  of  this  from 
the  scatter  plots  of  Figs. (16).  A  striking  feature  of  these  scatter  plots  is 
the  existence  of  a  large  region  on  the  underside  of  the  vortex,  which  is 
almost  evacuated  of  particles,  and  which  lies  along  the  vortex  separatrix. 
We  believe  that  the  evacuated  regions  in  Figs.(16)  are  formed  when  parti¬ 
cles  which  are  E  x  B  drifting  counter-clockwise  around  the  vortex  impact 
into  the  wall,  at  a  point  at  the  level  of  the  vortex  center.  The  result  is 
a  depletion  of  the  region  on  the  underside  of  the  vortex:  thus,  the  evac¬ 
uated  regions  can  be  considered  to  be  the  “wake”  inside  the  plasma,  of 
the  scrape-off  of  particles  right  at  the  wall.  Not  surprisingly,  the  existence 
of  the  evacuated  regions  indicates  that  particle  orbits  are  most  unstable 
on  the  separatrix,  a  generic  feature  of  dynamical  systems[39].  This  is  also 
qualitatively  confirmed  in  Figs. (23):  we  can  see  that  an  electron  moving 
right  on  the  separatrix  will  have,  at  intervals  of  a  bounce  period,  many 
opportunities  for  grazing  the  wall  (this  is  particularly  clear  for  Electron  3). 

7.2  Dynamics  without  Electron  Gyromotion 

The  results  of  the  previous  section,  in  which  we  saw  that  by  itself,  the  E  x  B 
motion  of  the  test  electrons  can  lead  to  diffusion  and  loss  of  confinement, 
suggested  that  we  run  a  simulation  in  which  the  gyro-motion  of  all  electrons 
was  suppressed,  so  as  to  demonstrate  that  only  the  E  x  B  motion  of  the 
electrons  is  important.  This  notion  is  certainly  reasonable,  in  view  of  the 
low-frequency  nature  of  the  spectra  shown  in  Figs. (20),  in  which  the  spectra 
were  confined  to  frequencies  |u;j  <C  u;ce. 

To  run  a  simulation  with  suppressed  gyromotion,  and  yet  maintain  a 
steady-state  loss  of  electrons  to  the  wall,  we  had  to  re-introduce  an  effect 
of  the  finite  gyroradius  by  creating  an  artificial  scrape-off  layer,  defined 
by  requiring  that  all  electrons  within  a  fixed  distance  x  <  pe  be  immedi¬ 
ately  absorbed  by  the  wall.  With  this  modification,  we  ran  a  simulation 
of  a  system  half  as  long  as  in  Run  1,  with  Ly  =  128  and  half  the  in¬ 
jection  rate,  but  with  otherwise  identical  parameters.  The  results  were 
qualitatively  very  similar  to  those  of  Rim  1,  with  the  system  establishing 
large-amplitude,  steady-state  vortices  (e  \8<f>\  /T,  ~  2),  and  with  a  transport 
which  maintained  the  particle  densities  at  values  very  close  to  those  in  Run 
1  (ne,j  «  2.1  peak  densities,  with  similar  density  profiles).  The  structure  of 


the  shear  layer  was  also  very  similar  to  that  found  in  Run  1. 

These  results  indicate  that  E  x  B  motion  dominates  the  electron  dy¬ 
namics,  and  suggest  that  future  simulation  studies  incorporate  an  electron 
mover  which  entirely  ignores  the  gyromotion  and  follows  only  the  E  x  B 
motion  of  the  guiding  centers.  This  modification  should  result  in  a  con¬ 
siderable  reduction  in  the  computation  times  (in  the  simulation  above,  we 
merely  damped  the  gyromotion,  and  were  still  limited  by  constraints  of  the 
explicit  time  step). 

7.3  Numerical  Estimates  of  the  Diffusion  Coefficient 

We  shall  now  estimate  the  diffusion  coefficient  using  several  different  mea¬ 
sures  of  transport.  Our  approach  is  qualitative,  and  aims  for  an  estimate 
of  the  average  diffusion  rate  across  the  system. 

From  Figs. (23),  we  can  estimate  a  characteristic,  “diffusive”  time 
for  the  diffusion  of  an  electron  inside  the  vortex.  We  define  Tdijf  as  the 
time  required  for  an  electron  to  migrate  to  the  separatrix.  Averaging  over 
the  five  test  electrons,  we  find: 

UciTdiif  =  215  (68) 

Note  that  this  is  the  time  for  the  diffusion  of  an  electron  across  the  vortex 
radius,  or  a  distance  roughly  Lxj 2.  We  can  now  compare  tj,//  to  a  global 
“confinement  time”  tj0„ ,  where  we  define  r;0„  to  be  the  average  time  any 
one  electron  spends  in  the  system.  Thus,  in  the  steady-state,  where  rate 
of  creation  equals  outward  flux,  we  must  have  s  =  Ne/Tiott,  where  s  is  the 
electron-ion  pair  creation  rate,  and  Ne  the  total  number  of  electrons.  With 
the  known  steady-state  values  Ne  =  16400  and  s  =  1.005,  we  immediately 
find: 

UciTio,.  =  410  (69) 

Thus,  we  find  that  r;0„  is  comparable  to  but  somewhat  greater,  with 

Tio,»/'Tdiff  =  1.9-  This  disparity  is  not  surprising,  insofar  as  r(o„  accounts 
for  the  migration  of  particles  across  the  entire  system,  while  rd,j j  was  an 
estimate  for  their  diffusion  across  roughly  half  the  distance  (and  only  in 
the  well  of  the  vortex). 


From  the  diffusion  time  given  in  Eq.(68),  we  can  in  turn  estimate  the 
diffusion  coefficient.  Since  Tdtjj  is  the  time  for  diffusion  across  a  distance 
Lx/2 ,  we  have  the  estimate: 


Di,if  = 


(LJ  2)2 


0.015 


where  the  superscript  indicates  an  estimate  based  on  test-particle  diffusion. 
How  does  the  value  given  above  compare  to  the  Bohm  diffusion  coefficient? 
The  latter  coefficient  is  given  by: 

rp  2 

Dx  =  a—  =  a —  (<1) 

eti  u;c> 

where  a  a;  1/16[40].  With  the  values  given  in  Table  1,.  we  find  that  Dx  a; 
0.0625,  so  that  Ddx,fi/D f  as  1/4. 

We  can  independently  estimate  the  diffusion  coefficient  by  requiring 
that  the  solution  of  the  diffusion  equation  be  compatible  with  the  observed 
density  profile.  The  diffusion  equation  is  given  by: 

d 

—n(x,t)  =  Dx—n(x,t)  +  <r  (72) 

wher^  a  =  s/LxLy  is  the  distributed  source  term  and  n(x,t)  the  y-averaged 
particle  density.  We  solve  Eq.(72)  in  the  steady-state  ( dn/di  =  0),  with 
boundary  conditions  n(x  =  0)  =  0  and  dn(x  =  Lx)/dx  =  0.  The  resulting 
solution  is  the  parabolic  profile: 


«(*)  =  Jd~x(2LX  ~  x)  (73) 

With  the  peak  density  known  from  the  simulation  ( n(Lx )  as  2.4  at  u ;c,t  = 
700,  see  Fig. (15)),  we  find  the  diffusion  coefficient: 


rjden*  ® ^ x 

1  “  2 n(Lx) 


0.026 


where  the  superscript  indicates  that  Dxens  is  derived  from  the  density  pro¬ 
file.  The  rough  agreement  between  the  estimates  of  Eqs.(70)  and  (74) 
indicates  that  our  calculations  are  consistent.  The  result  Dxen *  >  Dfff 
suggests  that  particles  injected  at  random  in  the  system  diffuse  outward 


faster  than  those  injected  right  in  the  center  of  the  vortex,  but  that  the 
disparity  in  diffusion  times  in  not  very  large. 


7.4  Estimate  of  Discrete-Particle  Effects 

We  shall  now  gauge  the  importance  of  collisional  effects  in  the  cross-field 
transport. 

Like  real  particles,  the  numerical  particles  undergo  scattering  and  diffu¬ 
sion  from  two  sources  of  discrete-particle  interactions:  short-range  binary 
collisions,  and  transport  by  long-wavelength,  thermal  fluctuations.  In  our 
range  of  parameters,  it  can  be  estimated  that  the  E  x  B  motion  induced 
by  the  thermal  background  is  the  dominant  contribution  to  the  electron 
diffusion[41]  (this  is  because  pe  <C  Xje).  The  resulting  two-dimensional 
diffusion  coefficient  is,  in  physical  units[42,41]: 

D‘,M  =  (l  +"jA4)'‘/’log(fcn,„i/2>r)  (75) 

In  Eq.(75),  the  numerical  value  of  the  constant  0  remains  uncertain.  The 
factor  3  is  a  constant  which  accounts  for  precise  definitions  of  correlation 
times,  and  for  numerical  smoothing  effects,  factors  which  are  not  accurately 
determined  in  the  analytic  derivation  of  Eq.(75).  The  essential  feature  of 
Eq.(75)  is,  however,  the  strong  dependence  of  the  collisional  diffusion  coef¬ 
ficient  on  n Xd.  If  the  transport  in  Run  1  is  mostly  collisional,  then  changing 
the  plasma  “graininess”  parameter  n\2d,  while  keeping  the  other  plasma  pa¬ 
rameters  fixed,  should  result  in  strong  changes  in  the  diffusion  coefficient, 
and  hence  in  markedly  different  steady-state  conditions:  thus,  we  have  a 
test  for  determining  whether  the  transport  predominantly  collisional. 

With  the  steady-state  parameters  of  Run  1,  upi/u jci  =  1.2,  TJ eB  =  1, 
and  a  numerically  imposed  kmax  =  0.237r/Ay,  we  find  the  estimate  for  the 
collisional  diffusion  coefficient: 


pc ou  2.24 3 

1  ^ 


(76) 


where  we  have  left  nX 2d  a  variable  parameter.  In  order  to  test  the  applicabil¬ 
ity  of  Eq.(76),  we  performed  a  series  of  simulations  in  which  we  varied  n\d, 
but  kept  uce  and  u>c,  fixed,  and  in  which  in  we  also  kept  the  same  creation 


of  charge  per  unit  time  from  one  simulation  to  the  next.  This  was  done  by 
changing  the  electric  charge  and  mass  per  numerical  particle  (which  were 
made  bigger  to  decrease  nXj),  while  keeping  e/me,  me/m,  and  the  rate  of 
injection  of  charge  e  x  s  =  (Charge  per  particle)  x  (Rate  of  pair  creation) 
fixed.  In  each  simulation,  we  determined  the  diffusion  coefficient  from 
Eq.(74): 


~Ml7) 

where  n(Lx )  is  the  central,  y-averaged  steady-state  density.  The  results  are 
shown  in  Fig. (27)  (Run  1  corresponds  to  \JnX2d  =  7.5).  In  Fig. (27),  it  can 
be  seen  that  the  observed  diffusion  coefficient  does  not  have  the  collisional 
dependence  implied  by  Eq.(76).  In  fact,  despite  strong  variations  in  nAy 
we  obtained  steady-states  with  very  similar  final  density  profiles  and  very 
similar  values  of  the  density-dependent  parameters  ay  and  ay.  Only  at  the 
smallest  value  of  n\2d,  nXd  =  4,  do  we  begin  to  see  indications  of  enhanced 
diffusion.  The  vortex  structures  also  remained  qualitatively  the  same,  but 
became  noisier  at  the  smaller  values  of  the  plasma  parameter.  These  results 
lead  us  to  believe  that  collisional  effects  are  not  the  dominant  transport 
mechanism  in  our  simulations. 


7.5  Analytic  Estimate  of  the  Diffusion  Coefficient 

In  the  previous  sections  we  have  argued  that  the  large  vortices  cannot 
account  uirectly  for  transport  in  the  system,  and  that  this  transport  does 
not  scale  as  would  be  expected  from  collisional  effects.  These  results  suggest 
that  the  turbulent  component  of  the  edge  spectrum  is  the  main  source  of 
transport,  and,  in  this  section,  we  shall  estimate  its  contribution  to  the 
diffusion  coefficient. 

We  shall  make  several  simplifying  assumptions.  We  assume  that  in  the 
steady-state  the  sheath  organizes  itself  so  as  to  have  a  thickness  /x,  and 
that  the  vortices  have  an  average  spacing  in  y  of  L,  with  diameter  lv.  We 


Figure  27:  Dependence  of  the  diffusion  coefficient  on  the  plasma  parameter 
nAj;  Run  1  has  (nA^)1^2  =  7.5. 


/„  =  c2Pi  (79) 

where  C\  and  c2  are  dimensionless  constants.  We  have  already  seen  that 
/„  ss  lx  «  5/3,-,  so  that  c2  ~  5,  and  we  expect  that  ly  will  be  somewhat  larger 
than  the  system  length  in  Run  1,  so  that  Cj  >  40.  We  also  assume  that 
each  vortex  has  an  average  depth  —<j>max,  such  that: 

,  _  .  T, 

tymax  —  ^*3 

e 

where  C3  is  another  dimensionless  constant,  with  C3  «  2.  In  what  follows,  we 
shall  assume  that  the  edge  turbulence  is  independent  of  the  edge  density,  at 
least  when  Upi  >  u;c,-.  This  assumption  is  motivated  by  the  observation  that 
the  cross-field  equations  (17,18)  appear  to  be  weakly  dependent  on  density 
when  either  Wp,  ujci  or  Wp,  <C  wc<,  because  in  both  limits  we  recover  the 
two-dimensional  Navier-Stokes  equation.  As  a  consequence,  we  take  the 
constants  Ci,  c2  and  C3,  as  well  as  other  parameters  characterizing  the  edge 
turbulence,  to  be  density-independent.  Note  also  that  by  scaling  all  lengths 
to  the  gyro-radius  as  in  Eqs.(78,79),  and  the  peak  amplitudes  to  the  ion 
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thermal  energy  as  in  Eq.(80),  we  have  essentially  stated  all  dependence  on 
magnetic  field  and  temperature  as  well. 

Let  us  now  consider  the  edge  spectrum  at  some  fixed  x,  say  x  =  lzj 2.  We 
shall  assume  that  the  spectrum  has  a  universal  dependence  on  wavenumber 
which  is  expressed  by: 

<  |^(L*y)|  >=  A(kypi)  <  |<£(0)|  >  (81) 

where  A  is  a  continuous  function,  with  >1(0)  =  1,  which  decreases  rapidly 
when  |Ly|p,  >  1.  In  Eq.(81),  we  have  the  Fourier  transform: 

4>(ky)  =  fV  0{y)e~ik^dy  (82) 

Jo 

In  particular,  we  have  the  estimate: 


< 


|^(0)|2  >=<  f  *  <f>{y)dy  WrU2  =  (— )  (  — )  l\  (83) 

J 0  \  C|  /  \  6  / 


Now  the  diffusion  coefficient  for  turbulent  x -directed  transport  is  given  by 
the  integral  over  the  velocity  autocorrelation  function: 

Dturb  =  <  vx(0 )vx(t)  >  dt  (84) 

Jo 

where  vx(t)  refers  to  the  velocity  of  a  test  particle.  At  a  given  (y,f),  the 
x-component  of  the  electron  velocity  is  given  by: 


vx  (y.t)  =  - 


B , 

B 


ly  ky 


(85) 


Introducing  (85)  into  Eq.(84),  we  can  then  decompose  the  diffusion  coeffi¬ 
cient  into  a  sum  over  the  contributions  from  separate  modes.  In  evaluating 
(84),  we  simplify  the  integration  over  time  by  introducing  a  correlation  time 
for  the  interaction  of  an  electron  with  a  given  wavenumber: 


rc(ky)  =  ucil  T(kypi)  (86) 

where  once  again  we  have  isolated  the  essential  physical  parameters,  and 
introduced  a  dimensionless  function  T(kyp, ),  such  that  T(l)  %  1,  and 


T(kyp, )  — *  0  as  kyp,  — »  oo.  We  also  introduce  a  low- wavenumber  cut¬ 
off  for  the  turbulent  spectrum,  by  excluding  from  the  contributions  to  the 
diffusion  coefficient  those  wavenumbers  with: 

kyPi  <  0  (87) 

where,  say,  0  ss  1/2.  Introducing  (85)  into  Eq.(84),  we  obtain: 

£‘ur6  =  f  E  kl  <  \<KK)\2  >  Tc(h)  (88) 

ly  kvPi>& 

Using  the  various  normalizations  (78,79,80),  we  obtain: 

Dtxurb  =  co  ^  (89) 

where  cq  is  the  dimensionless  constant  defined  by  the  expression: 

co  =( —V  -  r  dKK2A(K)T(K)  (90) 

\  C\  J  ir  J0 

Taking  cx  «  40,  c2  w  5  and  c3  «  2,  and  assuming  the  integral  to  be  of  order 
1,  we  find  an  estimate  for  c<j: 

c0  «  0.02  (91) 

a  value  somewhat  smaller  than  implied  by  Eq.(74),  but  of  the  right  order 
of  magnitude. 

Eq.(90),  which  predicts  a  Bohm-type  of  diffusion,  is  also  consistent  with 
the  general  result  of  Taylor  and  McNamara[43],  who  found  that  in  two- 
dimensions,  E  x  B  transport  will  always  lead  to  a  1  / B  scaling  of  the  diffusion 
coefficient.  Our  assumption  that  the  edge  turbulence  is  density-independent 
must  then  result  in  a  diffusion  coefficient  which,  dimensionally,  can  only 
have  a  Bohm-like  dependence  on  B  and  T,. 

For  the  specific  case  of  Run  1,  we  can  estimate  the  diffusion  coefficient 
from  the  discrete  sum  of  Eq.(84).  Assuming  that  0  =  0.5,  Tc(ky)  =  u>~1. 
and  using  the  numerical  values  of  the  spectra  as  obtained  in  Fig. (19),  we 
find  the  estimate: 

D‘xurb  =  0.04  (92) 


79 


a  very  rough  estimate,  which  nonetheless  is  in  accordance  with  the  numer¬ 
ical  value  obtained  in  Eq.(74).  We  shall  now  consider  the  density  depen¬ 
dence  of  the  diffusion  coefficient,  as  measured  in  our  numerical  simulations, 
thereby  qualifying  the  range  of  validity  of  Eq.(89),  which  assumes  density- 
independence. 


7.6  Scaling  of  the  Diffusion  Coefficient  with  Density 

We  shall  now  consider  the  effect  of  particle  density  on  the  diffusion  co¬ 
efficient,  as  measured  in  the  numerical  simulations.  We  performed  a  se¬ 
ries  of  numerical  simulations  which  were  variations  of  Run  1  (but  in  a 
smaller  system,  with  Ly  =  128),  over  which  we  studied  the  effect  of  vary¬ 
ing  the  final  Wp,/u>c,  ratio.  We  achieved  this  by  varying  the  creation  rate 
s  (the  rate  of  creation  of  electron-ion  pairs  per  time  step)  over  the  range 
0.157  <  s  <  5.02,  which  resulted  in  a  range  of  average  steady-state  densities 
such  that  0.63  <  <  8.07.  We  then  estimated  the  diffusion  coefficient 

using  Eq.(74): 

„  all 
1  2  n(Lx) 

where  a  =  s/ LxLy  is  the  injection  rate  per  unit  area.  The  results  are  shown 
in  Fig. (28),  where  we  plot  Dx  against  The  figure  indicates  that  at 

lower  densities,  in  the  range  u>^/u <  3,  the  diffusion  coefficient  has  an 
approximately  linear  dependence  on  the  density.  At  higher  densities,  with 
>  3,  the  diffusion  coefficient  levels  off,  and  assumes  what  appears 
to  be  a  density  independent  value  of: 


(93) 


Dx  =  0.0A^  (94) 

The  density-independence  of  the  diffusion  coefficient  is  consistent  with  the 
derivation  of  the  previous  section,  where  we  argued  that  provided  uipj 
1,  the  diffusion  coefficient  should  be  density  independent. 

In  the  simulations  leading  to  Fig. (28),  we  were  able  to  explore  only 
a  limited  range  of  densities,  amounting  to  an  order  of  magnitude.  This 
limitation  occured  because  we  desired  to  keep  <  1,  while  running 

the  simulation  with  the  small  numerical  mass  ratio  m,/me  =  40.  To  reach 
higher  steady-state  densities,  while  maintaining  <  1,  would  have 
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Figure  28:  Dependence  of  the  diffusion  coefficient  on  the  steady-state  den¬ 
sity. 

required  using  larger  mass  ratios,  and,  with  the  explicit  scheme  used  in 
ES2,  would  have  resulted  in  prohibitively  long  computer  runs.  As  noted  in 
Section  7.2,  these  limitations  should  be  overcome  in  a  future  modification 
of  ES2,  incorporating  an  electron  E  x  B  mover. 

Notwithstanding  the  limitations  in  the  range  of  densities  explored,  we 
expect  that  the  levelling-off  of  the  diffusion  coefficient  with  larger  den¬ 
sities  reflects  a  general  trend,  and  should  be  applicable  at  considerably 
higher  densities,  for  instance  in  a  fusion  environement,  where  one  might 
have  u-p, as  large  as  103. 

We  should  also  note  that  while  the  diffusion  coefficient  appears  to  be 
rather  insensitive  to  the  steady-state  density,  the  nature  of  the  edge  turbu¬ 
lence  does  appear  to  be  modified  by  the  larger  injection  rates.  In  Figs. (29), 
we  compare  potentials  for  two  different  injection  rates,  which  result  in 
steady-state  densities  such  that  iap,/u ;ct  —  2.1  and  u.’p,/u;c,  =  2.8.  The  qual¬ 
itative  difference  in  these  plots  suggests  that  for  uipt/u ;c,  >  2  there  occurs  a 
transition  to  a  more  turbulent  regime,  in  which  the  long-lived  vortices  seen 
in  Run  1.  wpi/u,'Ci  =  1.3.  are  replaced  by  less  stable  structures. 

Finally,  we  should  emphasize  that  the  diffusion  coefficient  of  Eq.(89)  is 
an  average  for  diffusion  in  the  sheath,  and  is  meant  to  be  nonzero  onlv  in  a 
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Figure  30:  Global  scenario  for  steady-state  turbulence  and  transport  in  the 
sheath. 

8  Summary 

In  Fig. (30),  we  have  summarized  the  interacting  mechanisms  which  ac¬ 
cording  to  our  scenario,  drive  the  turbulent,  steady-state  cross-field  sheath. 
The  scenario  proceeds  as  follows.  The  presence  of  the  wall  imposes  the 
formation  of  the  initial  edge  shear  layer,  which  is  unstable  to  the  Kelvin- 
Helmholtz  instability.  This  instability  is  in  turn  the  source  of  the  edge 
turbulence,  which  after  a  transient  phase,  reaches  a  steady-state.  In  this 
new,  dynamic  equilibrium,  the  long-wavelength  components  of  the  turbu¬ 
lence,  the  vortices,  modify  the  edge  conditions,  by  broadening  the  edge 
layer,  and  by  maintaining  large-amplitude,  drifting  potential  structures. 
Furthermore,  the  short-wavelength  components  of  the  turbulence  create  a 
cross-field  transport,  by  diffusing  the  electrons  out  of  the  system.  As  a 
result,  an  ambipolar,  outward  flux  of  electrons  and  ions  is  set  in  motion, 
with  a  diffusion  coefficient  which  scales  as  the  Bohm  diffusion  coefficient. 


Thus,  with  a  constant  injection  of  electron-ion  pairs,  the  system  maintains 
a  steady-state,  in  which  the  linear  edge  instability,  the  nonlinear  dynamics 
of  the  turbulence,  and  the  outward  particle  diffusion  all  balance  each  other. 

Our  principal  results  can  be  summarized  as  follows: 

1.  The  y- uniform  cross-field  sheath  is  linearly  unstable  to  the  Kelvin- 
Helmholtz  instability,  with  growth  rates  predicted  by  fluid  theory  in 
the  range  \ky\pi  <C  1. 

2.  The  turbulent  cross-field  sheath  has  an  equilibrium  thickness  of  order 
lx  ~  5 Pi. 

3.  The  edge  turbulence  principally  consists  of  large- amplitude  vortices, 
with  \S(j>\  ~  3T,/e,  and  with  a  spacing  of  at  least  ly  ~  40 p,. 

4.  The  drift  velocity  of  the  vortices  parallel  to  the  wall  is  of  order  vt,/2, 
leading  to  a  turbulent  spectrum  with  dispersion  relation  at  small 
wavenumbers  |u;|  \ky\vti/2. 

5.  The  edge  spectrum  extends  in  wavenumber  to  |k[p,  «  1,  beyond 
which  it  rapidly  decreases;  the  shorter-wavelength  components  have 
a  broad  spectrum  in  frequency,  with  maximum  width  |<5u>|  ~  u>ci. 

6.  Transport  in  the  sheath  is  Bohm-like,  and  is  consistent  with  a  model 
in  which  the  electrons  are  driven-out  by  the  short-wavelength  turbu¬ 
lence,  with  a  subsequent,  ambipolar  diffusion  of  the  ions.  An  estimate 
of  the  diffusion  coefficient  is: 

Dxrb  =  0-04^,  ujpi  >  2uc,  (97) 

At  lower  edge  densities,  such  that  wp;  <  2u ;c,-,  the  expression  for  the 
diffusion  coefficient  has  an  additional  factor,  which  is  approximately 
proportional  to  the  density. 

7.  For  a  given  edge  density  n0,  the  sheath  induces  a  particle  flux  to  the 
wall  which  is  of  order: 

rx  =  10  2u<,Uo,  iopi  >  2u>cj  (08) 
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Finally,  let  us  say  a  few  words  regarding  similarities  and  differences 
between  our  model  and  that  of  Horton  et  al. [16] .  In  [16],  an  implicit  parti¬ 
cle  code  is  used  to  explore  the  evolution  of  an  initially  symmetric  velocity 
shear  layer.  The  evolution  of  the  Kelvin-Helmholtz  instability  in  this  con¬ 
figuration  is  similar  to  that  which  we  found  for  the  cross-field  sheath:  the 
instability  grows,  saturates  into  large-amplitude  vortices,  which  then  un¬ 
dergo  multiple  coalescence,  and  finally  merge  into  one  large  vortex.  There 
are  nonetheless  basic  differences  between  the  model  of  [16]  and  our  model: 

1.  In  our  model,  the  primary  concern  is  plasma- wall  interaction;  in  ES2 
the  boundary  is  fully  absorbing,  while  the  boundaries  in  [16]  are 
purely  reflecting,  and  are  placed  far  from  the  shear  layer. 

2.  Because  of  the  presence  of  the  wall,  the  cross-field  sheath  and  its 
vortices  form  in  a  completely  spontaneous  and  self-consistent  manner. 
On  the  other  hand,  the  model  of  [16]  requires  an  external  electric  field 
in  order  to  initiate  the  shear  layer. 

3.  There  is  a  basic  difference  in  space  and  time  scales.  Horton  et  al. 
consider  a  shear  layer  with  width  a  »  p,,  and  their  imphcit  particle 
code  is  intrinsically  limited  to  low  frequencies,  u ;  <C  u>CJ.  On  the 
other  hand,  our  model  includes  much  shorter  space  and  time  scales: 
the  sheath  develops  a  thickness  of  a  few  gyro-radii  across,  and  the 
turbulence  extends  to  frequencies  with  w  »  wci. 

4.  In  our  model,  the  cross- field  sheath  exhibits  strong  density  variations, 
whereas  [16]  considers  a  shear  layer  with  a  constant  density. 

Thus,  despite  overlap  with  our  model,  especially  as  regards  the  fluid 
evolution  of  the  Kelvin-Helmholtz  instability,  the  model  of  [16]  nonetheless 
applies  to  a  very  different  range  of  physical  parameters. 


9  Conclusions 

Our  numerical  simulations  have  shown  that  the  cross-field  sheath  between  a 
wall  and  a  plasma  is  not  a  static  structure,  but  is  in  fact  a  turbulent  bound¬ 
ary  layer,  with  strong  potential  fluctuations  and  anomalous  particle  trans¬ 
port.  The  driving  mechanism  for  this  turbulence  is  the  Kelvin-Helmholtz 
instability  which  arises  from  the  sheared  particle  drifts  created  near  the  wall 
by  the  strongly  noneutral  sheath.  Provided  it  is  replenished  by  an  internal 
flux  of  particles,  coming  for  instance  from  a  central  bulk  plasma  or  from  a 
diffuse  ionization  of  neutrals,  the  sheath  will  maintain  itself  in  a  dynamic 
equilibrium,  in  which  the  linear  edge  instability,  the  nonlinear  dynamics 
of  the  particles  and  the  outward  particle  diffusion  all  balance  each  other. 
It  is  important  to  emphasize  that  the  turbulent  behavior  of  the  sheath  is 
a  completely  spontaneous  phenomenon,  which  arises  from  the  plasma-wall 
interaction,  and  which  does  not  require  the  imposition  of  external  fields  for 
its  sustenance. 

We  have  found  that  the  cross-field  sheath  assumes  an  equilibrium  thick¬ 
ness  of  order  lx  ~  5  p,,  and  that  it  maintains  large  vortices,  with  amplitudes 
6(f)  ~  T,/e,  which  drift  parallel  to  the  wall  at  roughly  half  the  ion  thermal 
velocity.  The  sheath  also  maintains  a  large,  spatially-averaged  potential 
drop  from  the  wall  to  the  plasma,  with  A4>  «  —1.5 T,/e,  in  sharp  distinc¬ 
tion  with  the  unmagnetized  sheath,  where  the  plasma  potential  is  higher 
than  at  the  wall.  The  average  velocity  shear  profile  remains  linearly  unsta¬ 
ble  at  all  times,  a  feature  which  shows  that  the  saturation  mechanism  is 
not  quasilinear  (relaxation  of  the  space-averaged  shear  profile),  but  strongly 
nonlinear. 

A  central  result  of  this  paper  is  that  the  sheath  induces  an  anomalous 
transport  of  particles,  which  is  due  to  short-wavelength  turbulence,  and 
which  scales  very  much  like  Bohm  diffusion,  at  least  when  ujpi  >  2 uc{.  At 
lower  densities,  such  that  Wp,  <  2 the  diffusion  coefficient  is  fotrnd  to 
have  an  additional  factor,  proportional  to  the  density.  These  results  enable 
us  to  model  the  entire  cross-field  sheath  by  a  simple  boundary  condition, 
which  relates  the  particle  flux  through  the  sheath  to  the  edge  density.  This 
boundary  condition,  which  simply  measures  the  sheath’s  “impedance”  to 
particle  flow,  should  be  very  valuable  as  input  in  any  model  designed  to  ob¬ 
tain  the  bulk  plasma  properties,  and  in  which  the  detailed  sheath  dynamics 
are  unimportant. 
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We  believe  that  future  numerical  simulations  of  the  cross-field  sheath 
should  explore  a  larger  range  of  parameters,  in  system  size  and  in  density, 
than  were  feasible  in  the  present  work,  in  order  to  verify  the  scaling  laws 
that  we  found  to  be  valid  for  a  more  limited  range  of  parameters.  For 
instance,  we  were  limited  to  systems  in  which  a  single  vortex  survived  in 
the  steady-state,  while  we  believe  that  longer  systems  should  be  able  to 
accomodate  many  vortices.  Similarly,  we  were  limited  to  exploring  a  range 
of  plasma  densities  such  that  <  10,  while  we  would  like  to  explore 

values  as  high  as  u^/u =  103,  representative  of  a  fusion  environment.  A 
relatively  straightforward  modification  of  our  simulation  code  ES2,  to  allow 
for  an  E  x  B  electron  mover,  should  make  all  of  these  extended  studies 
feasible  at  a  reasonable  cost  in  computer  time. 

Other  avenues  of  research  on  cross-field  dynamics  might  involve  study¬ 
ing  the  behavior  of  a  free  plasma  boundary,  far  from  the  wall,  in  order  to 
observe  the  relaxation  of  the  edge  layer.  A  central  question  here  is  whether 
the  Kelvin-Helmholtz  turbulence  can  induce  transport  over  very  many  ion 
gyro-radii,  or  whether  it  is  limited  to  a  “natural”  boundary  layer,  of  thick¬ 
ness  comparable  to  the  plasma  wall  sheath,  that  is  roughly  lx  ~  5 pt.  Let  us 
also  note  that  in  connection  with  active  plasma  devices,  one  might  study 
the  behavior  of  the  cross-field  sheath  when  an  external  electrostatic  po¬ 
tential  is  applied:  the  resulting  dynamics  should  bear  much  ressemblance 
to  those  of  the  one-species  magnetron,  but  with  the  distinct  physics  of  a 
two-species  plasma.  Finally,  we  note  that  our  simulation  results  also  raise 
more  fundamental  questions  regarding  the  nature  of  the  edge  turbulence, 
which  was  only  partially  explained  by  the  fluid  theory.  This  suggests  fur¬ 
ther  analytical  research,  applying  a  fully  kinetic  treatment  to  the  problem, 
with  perhaps  specially  designed  simulations  to  verify,  for  instance,  vortex 
stability  and  the  mechanisms  for  generating  short-wavelength  fluctuations. 


Acknowledgments 

The  author  wishes  to  thank  Prof.C.K.Birdsall  for  his  contributions  to  this 
research,  through  various  discussions,  and  especially  for  initiating  this  work, 
by  suggesting  the  possibility  of  the  Kelvin- Helmholtz  instability  in  the  two- 
dimensional,  cross-field  sheath.  The  author  also  wishes  to  thank  Dr.Liu 
Chen  for  many  detailed  and  stimulating  exchanges,  and  Dr.M.J.Gerver  and 
Prof.A.J.Lichtenberg  for  various  illuminating  discussions.  Finally,  thanks 
go  to  Mr. 0. Garcia  of  the  Electronics  Research  Laboratory,  for  his  assistance 
with  the  graphical  work.  This  research  was  supported  by  U.S.  Department 
of  Energy  Contract  No.FG03-86ER53220,  by  U.S.  Office  of  Naval  Research 
Contract  No.  N14-80-C-0507  and  by  a  MICRO  grant  with  a  gift  from  the 
Varian  Corporation. 


Appendix 

A  The  One-Dimensional  Sheath 


In  this  appendix,  we  study  the  structure  of  the  sheath  when  the  geometry 
is  restricted  to  a  single  spatial  dimension,  the  perpendicular  distance  from 
the  wall.  This  analysis  is  conceived  as  an  initial  value  problem:  we  consider 
the  sheath  that  forms  from  the  scrape-off  of  particles  in  an  initially  uniform 
plasma,  which  is  put  into  contact  with  the  wall  at  some  t  =  0.  The  result, 
for  t  2tt/u is  a  static  sheath  which  extends  some  distance  p,  into  the 
plasma.  In  our  discussion  of  this  initial- value  problem,  we  omit  treating  the 
effects  of  collisionality  and  of  a  distributed  source  of  particles,  even  though 
it  is  only  in  the  presence  of  these  processes  that  the  one- dimensional  sheath 
can  maintain  a  steady-state  transport  to  the  wall.  Such  an  analysis  has  been 
undertaken  by  several  authors[44,45].  In  the  present  Section,  our  goal  is 
less  ambitious.  We  consider  the  formation  of  the  one-dimensional  sheath 
as  a  transient  prelude  to  two-dimensional  behavior,  and  as  such  restrict 
our  analysis  to  the  initial  sheath  formation.  The  results  of  our  study  make 
analytically  plausible  the  values  of  y-uniform  edge  potential  drop  and  of 
the  shear  layer  which  almost  instantaneously  appear  in  the  two-dimensional 
simulations. 

A.l  Analytic  Expressions  for  the  Electric  Field 

We  should  stress  that  the  analysis  in  this  Section  is  semi-quantitative.  This 
is  because,  despite  the  apparent  simplicity  of  the  formulation  of  the  initial- 
value  problem,  we  have  not  been  able  to  obtain  a  fully  self-consistent  so¬ 
lution  for  the  formation  and  structure  of  the  cross-field  sheath.  This  dif¬ 
ficulty  is  Unked  to  the  multi-dimensional  nature  of  the  problem,  where  we 
have  to  solve  for  the  time-evolution  of  a  distribution  function  f(x ,  t>r,  vy,  t) 
which  is  three-dimensional  in  phase  space.  Furthermore,  the  wall  boundary 
condition  is  such  that  if  /  is  Maxwellian  far  from  the  wall,  it  is  nonethe¬ 
less  strongly  non-Maxwellian  in  the  edge  layer,  and  thus  a  perturbation 
approach  is  not  feasible.  Another  complication  arises  because  (he  sheath 
thickness  is  of  the  same  order  as  the  gyroradius,  and  this  implies  that  /  has 
a  strong  dependence  on  all  of  its  independent  variables  (  /  varies  a  length 


scale  of  order  p,,  over  a  velocity  scale  vti  and  over  a  time-scale  u>~1 ).  In  short, 
there  is  no  small  parameter  to  facilitate  an  analytic  solution,  except,  as  we 
shall  see,  for  the  ratio  c u^/u^,  which  can  be  considered  small  in  a  fusion 
environment.  Paradoxically,  some  aspects  of  the  two-dimensional  behavior 
of  the  sheath  are  more  amenable  to  analysis  than  in  the  one-dimensional 
configuration. 

In  what  follows,  we  shall  assume  that  the  electron  gyro-radius  is  negli¬ 
gible  and  take  pe  — *  0.  Thus,  the  electrons  axe  “frozen”  onto  the  magnetic 
field  lines,  and  their  density  remains  constant  in  space  and  time.  With  this 
approximation,  we  need  to  consider  only  the  evolution  of  the  ion  distribu¬ 
tion  function,  f,(x,  vx,  vy ,  t),  in  the  presence  of  a  self-consistent  electric  field 
Ex(x,t),  taking  into  account  the  loss  of  ions  by  impact  to  the  wall.  Fur¬ 
thermore,  we  shall  consider  only  the  time-asymptotic  state  of  the  sheath, 
when  t  >■  2n juc{. 

With  the  wall  at  x  =  0,  we  assume  an  initially  uniform  and  Maxwellian 
ion  distribution  function: 

f,(x,vx,v tf,f  =  0)  =  ^^e_(u*+u»)/2u'S  0  <  x  <  co,  (99) 

To  bring  out  the  essential  dependencies  in  the  problem,  we  normalize  all 
variables  according  to: 


t'  =  ucit  x'  —  x/ pi,  (100) 

/'  =  v?J,  E'  =  Ex/B0vti ,  (101) 


where  p,  =  vt,/u is  the  ion  gyro-radius.  With  these  normalizations, 
the  Vlasov- Poisson  system  can  be  written  (dropping  the  primes  for  con¬ 
venience) 


d± 

dt 


dx 


2L 

dvT 


~§x  ~  ®  (/  ^ X ’  Vx'  Vy'  dVx  dVy  ~  O  ’ 


df 

(102) 

-1)' 

(103) 

with  the  initial  condition,  f(x,vx,vy,  0)  =  exp[— (v\  +  vy)/2]/2n.  The  sys¬ 
tem  of  Eqs.(102,103)  is  to  be  solved  for  0  <  x  <  oo.  In  Eq.(103),  we 
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have  defined  Q  =  ui^ / u*t  where  upi  is  the  plasma  frequency  of  the  uniform 
plasma  far  from  the  wall.  The  parameter  Q  is  the  only  remaining  indepen¬ 
dent  parameter  in  the  problem,  which  now  consists  in  finding  f(x,  vx,  vy ,  t) 
and  E(x,t)  for  t  — +  oo,  as  a  function  of  Q.  For  fusion  applications,  we 
have  Q>1,  and  this  will  afford  some  simplifications  in  our  analysis  of  the 
sheath  structure. 

Despite  their  apparent  simplicity,  Eqs.(102,103)  have  no  obvious  ana¬ 
lytic  solutions.  This  has  led  us  to  resort  to  a  semi-qualitative  approach 
is  solving  them.  Our  basic  assumption  is  to  make  what  might  be  termed 
a  “local”  approximation,  in  which  the  ion  density  at  any  one  point  is  an 
explicit  function  of  the  electric  field  evaluated  at  the  same  point.  We  first 
assume  that  at  any  distance  x,  the  ion  density  can  be  calculated  by  assum¬ 
ing  that  the  electric  field  at  that  point  has  been  constant  in  both  space 
and  time  since  t  =  0.  With  a  constant  electric  field,  an  ion  initially  at  a 
distance  x  and  with  velocity  (vx,vy)  will  not  collide  with  the  wall  provided 
that  its  velocity  lies  in  the  region: 

vy  >  -E  +  ^-(vl  -  x2)  =  g(x,vx),  (104) 

We  shall  now  approximate  the  local  “free”  ion  density  by  an  integral  which 
removes  all  particles  which  are  not  in  the  region  specified  by  (104): 

nF(E,x)  =  /  f\t(x,vx,vy)  dvz  dvy,  (105) 

Jvv>g(i,vx) 

Eq.(105)  does  not  embody  consistent  approximations,  and  in  fact  as¬ 
sumes  two  opposite  physical  limits,  one  for  the  ions  which  have  impacted 
into  the  wall,  and  another  for  the  “free”  ions  which  have  remained  in  the 
plasma.  For  the  unconfined  ions,  finite  orbit  effects  are  implicitly  taken 
into  account  by  assuming  that  they  have  indeed  impacted  into  the  wall, 
and  thus  must  be  removed  from  the  integration  in  Eq.(105).  On  the  other 
hand,  we  have  neglected  finite  orbit  effects  for  the  confined  ions,  because  in 
Eq.(105)  we  take  the  integrand  to  be  the  initial  Maxwellian,  without  con¬ 
sidering  the  mixing  of  contributions  from  different  areas  of  phase  space,  as 
would  be  required  in  a  calculation  using  the  full  ion  orbits.  Despite  these 
simplifications,  we  believe  that  the  spirit  of  the  calculation  is  correct,  in 
that  it  should  bring  out  the  proper  space  scales  and  field  amplitudes. 

Using  Eq.(109),  one  obtains  an  expression  for  the  local  free  ion  density: 
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nF  =  (2? lyf*  Le-x/2  dVy  e  U“/2erf  +  +  x2/2^1/2]  (106) 

where  “erf”  denotes  the  error  function.  For  the  analysis  that  follows, 
Eq.(106)  is  a  rather  untractable  expression  for  n/r,  and,  in  the  spirit  of 
the  previous  approximations,  we  shall  greatly  simplify  this  equation.  We 
first  make  the  hypothesis  that  the  fields  in  the  sheath  are  such  that  we  can 
everywhere  assume  that  either  E  1  (near  the  wall,  with  x  <  1),  or  that 
x  1  (in  the  tail  of  the  sheath,  where  E  <  1).  We  deliberately  neglect 
the  transition  region,  where  E  ~  1  and  x  ~  1,  hoping  that  the  essential 
behavior  comes  from  the  asymptotic  regions  near  and  far  from  the  wall.  As 
we  shall  see,  these  approximations  are  plausible  provided  Q  =  Wp,/u>^  ^  1. 

With  the  simplifying  assumptions  of  the  previous  paragraph,  we  can 
uniformly  neglect  the  ^-dependence  in  the  argument  for  the  error  function, 
and  furthermore,  we  set  the  lower  limit  of  the  integral  in  Eq.(lOC)  to  — oo. 
The  result  of  these  simplifications  is: 

np  —  erf  {(xE  +  x2/2)^2]  ( 107) 

We  now  impose  another  sweeping  approximation,  by  replacing  the  error 
function  in  Eq.(107)  by  the  expression: 


erf(z)  ~  1  —  e  z* ,  (108) 

an  approximation  which  incorporates  the  vanishing  of  the  error  function  at 
z  =  0  (albeit  with  the  wrong  dependence  on  z),  and  the  correct  exponential 
dependence  for  z  — >  oo. 

Given  np(E,x)  as  defined  by  Eq.(105),  we  then  introduce  it  into  Pois¬ 
son’s  equation: 

dE 

—  =  Q(nF-l),  (109) 

Using  (108)  in  Eq.(107),  and  introducing  the  resulting  np  into  Poisson's 
equation,  we  obtain: 


—  =  Q  exp  (-{xE  +  x2/2)) 
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(HO) 


This  equation  is  to  be  solved  with  the  sole  boundary  condition  that  E(x) 
vanishes  for  large  x,  a  condition  which  uniquely  determines  the  solution 
with  Q  as  a  parameter.  In  particular,  Eq.(llO)  uniquely  determines  the 
value  of  the  electric  field  at  the  wall,  E( 0)  as  a  function  of  Q. 

Despite  all  our  manipulations,  Eq.(16)  is  not  yet  analytically  solvable. 
However,  one  can  extract  important  information  on  the  nature  of  the  sheath 
by  examining  the  scaling  of  the  variables  on  Q.  Let  us  write  £  =  Ql?2x  and 
E  =  E/Ql/2.  We  the  have: 

^7  =  exp  (-(££  +  i2 /2Q])  (111) 

We  now  see  that  if  we  have  Q  1,  there  exists  an  edge  layer  extending  to 
some  £  ~  1  where  Eq.(lll)  is  independent  of  Q ,  and  has  the  form; 

^  =  exp  (~£E)  (112) 

In  physical  units,  this  edge  region  extends  to  x  &  p,/Q1^2  =  Adl,  and  the 
electric  field  has  an  amplitude  Ex  ~  Ql/2vt,  =  (vpi/wCi)vti  >  vtt.  With 
the  requirement  that  the  electric  field  vanishes  far  from  the  wall,  we  can 
roughly  estimate  E( 0)  from  Eq.(112)  by  assuming  that  in  the  right-hand 
side  of  the  equation,  we  can  take  the  electric  field  to  be  constant  and  equal 
to  its  value  right  at  the  wall.  Integrating  from  0  to  oo,  we  then  have  the 
identity: 

~  r°° 

E(0)=  e-^(0><f£  (113) 

Jo 

which  immediately  yields  E( 0)  =  1.  In  physical  units: 

£*(0)  ~  B0vtlH,  (114) 

U>CJ 

The  solution  of  Eq.(lll)  must  be  connected  to  the  “outer  solution"  of 
Eq.(112),  in  the  region  ^  >  1.  The  scaling  of  Eq.(lll)  breaks-down  when 
E(£)  becomes  so  small  that  it  becomes  comparable  to  the  second  term  in 
the  exponential  in  Eq.(lll).  This  occurs  when  x  ~  1  and  E  ~  1  /Q,  or  in 
physical  units,  we  have  Ex( Aj,)  ~  BqVh.  Beyond  this  point,  Eq.(16)  predicts 
the  decay  of  the  electric  over  a  length  comparable  to  the  gyro-radius  p,. 


To  summarize  the  results  of  our  approximations,  we  have  found  that 
with  u'px /u.'c,  >  1,  Poisson’s  equation  predicts  the  existence  of  a  “composite” 
sheath,  with  an  inner  and  an  outer  layer.  There  is  a  first  region,  an  “inner 
sheath",  extending  up  to  a  Debye  length  A*,  for  which  the  peak  electric 
field  was  given  by  Eq.(114),  Ex  ~  BoVt,(u)pi/iJct).  Beyond  this  inner  sheath, 
there  is  a  boundary  layer  extending  to  roughly  a  gyro-radius  p,  from  the 
wall,  with  a  peak  electric  field  of  order  Ex  ~  Bovtt. 

A. 2  Comparison  with  one-dimensional  particle  simu¬ 

lations 

We  now  briefly  compare  the  analytic  results  of  the  previous  Section  with 
those  of  one-dimensional  particle  simulations.  In  Fig. (31a, b),  we  show  the 
electric  fields  obtained  one-dimensional  particle-in-cell  simulations  of  the 
cross-field  sheath,  for  u>p,/uc,  =  1.9  and  6.3.  In  effect,  the  simulation  algo¬ 
rithm  exactly  solves  Eq.(102),  to  within  the  bounds  imposed  by  the  thermal 
noise  induced  by  the  finite  particles.  In  Figs. (31a, b),  we  can  see  that  the 
strict  division  of  the  sheath  into  inner  and  outer  layers,  as  was  implied  by 
the  considerations  of  the  previous  subsection,  is  not  observed  in  the  simu¬ 
lation  electric  fields  Ex(x).  On  the  other  hand,  the  overall  sheath  is  indeed 
an  ion  gyro-radius  in  thickness.  Furthermore,  we  have  found  that  the  scal¬ 
ing  for  the  peak  electric  field  at  the  wall  does  obey  Eq.(??)  remarquably 
well,  whenever  we  have  ujp,/uc,  >  1.  This  can  be  seen  in  Fig. (32).  where  we 
plot  the  value  of  the  electric  field  at  the  wall  as  a  function  of  the  parameter 
upi/u)ci.  Note  that  for  u>Pi/ojc,  <  1,  the  electric  field  appears  to  scale  as 
Ex{0)/Bovt,  ~  u£/u£. 
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Figure  31:  Profiles  of  the  cross-field  sheath  in  one-dimensional  simulations. 


Figure  32:  Scaling  with  density  of  the  electric  field  at  the  wall  in  a  one¬ 
dimensional  simulation  (full  dots).  Line  A  denotes  the  dependence  ujp,/uct< 
line  B  the  dependence 
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